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1. Introduction. In a recent paper’ Alexander Weinstein reconsidered axially sym- 
metric potential flows of an ideal, incompressible fluid and in particular discussed the 
flow generated by various axially symmetric distributions of sources. Weinstein based 
his analysis on the classical’ representation of the corresponding velocity potentials and 
stream functions in terms of improper integrals involving Bessel functions. In the present 
paper, which was suggested by and is complementary to Weinstein’s work, the potentials, 
stream functions, and velocity components for basic axially symmetric distributions of 
sources or vorticity over the circumference or the interior of a circle, are established in 
terms of elliptic integrals of the first and second kind. 

The relative merit of the alternative representation used here appears to be twofold. 
First, the approach via elliptic integrals yields an analytically more transparent de- 
scription than that afforded by the representation through discontinuous integrals of 
Bessel functions. In particular, further insight is gained into the cyclic character’ of 
Stoke’s stream function or the velocity potential in cases where the distribution of singu- 
larities is such as to give rise to a multiply connected domain of regularity in the me- 
ridional half-plane. The results here referred to are of course in complete agreement 
with those of Weinstein who originally clarified this aspect in connection with the stream 
function for the source ring. Secondly, the use of functions which have been tabulated 
exhaustively facilitates the physical interpretation of the results and renders possible 
the complete determination of the corresponding streamline patterns. By superposition 
of the foregoing basic axially symmetric flows and appropriately chosen uniform streams 
one may obtain a variety of technically significant flows around solid and annular shaped 
bodies and half-bodies of revolution. A comprehensive study of the various body shapes 
and associated streamline patterns so obtainable is currently in progress in cooperation 
with V. L. Streeter and P. C. Chu, of Illinois Institute of Technology, who have inde- 
pendently completed by aid of numerical integrations several flow patterns based on 
the homogeneous source disc. 

It should be recalled at this place that, historically, elliptic integrals were introduced 
early in connection with problems of rotationa] symmetry in potential theory. Thus 
the simple formula for the velocity potential of a homogeneous source ring appears in 
classical treatises as the gravitational potential of a homogeneous circumference, and 
the electro-magnetic analogue of the stream function for the vortex ring was given by 
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Maxwell® in terms of the complete integrals of the first and second kind. It appears, 
however, that the complete integral of the third kind and its relation to the incomplete 
integral of the second kind and Jacobian elliptic functions, have been neglected in the 
problems under consideration. In this connection it is interesting to note that G. M. 
Minchin* obtained an expression involving integrals of the third kind for what corre- 
sponds to the velocity potential of the vortex ring, but failing to apply the foregoing 
reductions, gives preference to the open representation in terms of a series of spherical 
harmonics. 

Finally, the authors understand that A. Van Tuyl, of the Naval Ordinance Labo- 
ratory, by formal transformation of the corresponding Bessel integrals has recently 
reached results which partly overlap with those given here. 

2. The governing equations. For the sake of convenience we cite here the basic 
equations governing the steady irrotational flow of an ideal incompressible fluid in the 
presence of axial symmetry. Referring the motion to the cylindrical coordinates (x, p, 6) 
where the x-axis is assumed coincident with the axis of symmetry, there exists a velocity 
potential (2, p) which satisfies Laplace’s equation 

0 0¢\ , @ Oo | : : 

ar (p =) + ap (p at = G (2.1) 
at all non-singular points of the field. Equation (2.1) may be considered as the condition 
of integrability assuring the existence of Stokes’ stream function y¥(z, p) as defined by 


Ox “a p’ Op Ox 


Choosing the additive constant in ¥(z, p) such that 


we can represent ¥(x, p) by the line integral 


e(z,p F," 4 ; 
¥(z, p) = | | -» 6 dx + » % ap | (2.3) 


¥(-2,0 


The choice of the path of integration is immaterial in simply connected fields. 


9 ] 4] 
= (3 ay += (2 ov) = 0 (2.4) 
Ox \p Ox Op \p Op 


The axial and radial velocity components are then given by 


From (2.2) follows 


t 


dp = l Oy — dp - ony 1 oy (2 5) 


. == sar : : 
‘ Ox p Op Op p Ox 


The curves ¥(z, p) = constant are the streamlines in the meridional half-plane and 


2r(¥. — v;) constitutes the total flow between the streamsurfaces generated by y = y, , 
y=". 
As was emphasized by Weinstein and is apparent from the definition (2.3), a given 
3J. C. Maxwell, Electricity and magnetism, Clarendon Press, Oxford, 1892, p. 339. 
‘G. M. Minchin, The magnetic field of a circular current, Phil. Mag., 35, 354-365, 1893. 
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velocity potential ¢(z, p) whose derivatives are single-valued guarantees the single- 
valuedness of the associated stream function ¥(z, p) if the domain of regularity in the 
half-plane p > 0 is simply connected. This condition, however, is in general not satisfied 
for source distributions off the axis of symmetry. 


p 


= 







iso-modular circle 
k=const. 


P(x,9) 








Ze=k' Z¢:-x' 








8'(0,-b) 
Fia. 1: Position Parameters 
3. Velocity potential and velocity components of the source ring. Consider a homo- 


geneous distribution of sources of total strength m (total flux 4rm) along the circle 
xz = 0, p = b. The distance R from a point P(z, p, 0) of the field to a point Q(0, b, 6) on 


the source ring is 





R = (2? + BD? + p’ — 2bp cos 8)” (3.1) 
= 7r,(1 — k’ sin’ g)'” (3.2) 
in which 
x— 0 T T 
va oy 9 (3.4) 
ri ri 


ri = 2° +(p + DB)’ 
(3.5) 


r; = x + (p — by’ 
so that r, , r. denote the distances B’P and BP respectively in Fig. 1. 
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The velocity potential at P is given by 


m [°" dé ; 
(x, p) = -= | : (3.6) 
? 2r /0 R 
To rationalize R, we introduce u by means of 
sin g = sn (u, k) 
—K<u<K (3.7) 
cos g = CN WU, k) 
which gives 
dy du 1 dé 
R =r, dn (u, k), — = =— (3.8) 
R i, 2R 
and thus 
»K 9 4 
m 2mKk ; 
o(z,p) = —-— | du=— (3.9) 
Tr) J-K TT; 


where K denotes the complete elliptic integral of the first kind for the modulus k. Pro- 


ceeding to the limit as k — 0, we confirm that along the z-axis 
= om 
¢(a, 0) = — (3.10) 
To 
where, 
ff, = r,(atp = QJ) = r,(atp = 0) = (2° + 2 iia (3.11) 
The components of the velocity field generated by the potential (3.9) are obtained 
from (2.5) by aid of the identity” 
dk l pre 
—- = —~73 [E — (k’)” K (3.12 
dk k(k’)° ! ° ) 
Here k’ is the complementary modulus, 
/ 2) 1/2 le 215 
k’ = (1 — k’)”* =- (3.13) 
rT; 
and E stands for the complete elliptic integral of the second kind referred to the modulus 
k. The computation yields 
2mz ~, 
ae E 


vz, p) = 3 
mT 2 
(3.14) 
. —0., 
v(x, p) = a (x = £5 — E) | 
) 
and, in particular, 

v(x, 0) = se v,(z, 0) = 0 (3.15) 


0 


‘For this and subsequent identities associated with elliptic functions and integrals, see for example 
E. T. Whittaker and G. N. Watson, Modern analysis, 4th ed., University Press, Cambridge, 1935, Ch. 


XXII. 
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The behavior of ¢, v, , v, in the vicinity of the source ring, i.e., in the neighborhood of 
the singular point B(0, b), is characterized by the following limits as r, and hence k’ 
tend to zero: 


m sin ¥ 


6+ log i 0, v. + — 0, ime = +6. (3.16) 
ab fs tbr. 


br. 
Here y denotes the polar angle at B (Fig. 1). Thus, 


; —2x —l 
sin y = a, cosy = & ~ -, —rsysr (3.17) 
2 2 





4. Stream function of the source ring. The stream function associated with the 
velocity potential (3.9) may be determined directly by virtue of its intrinsic hydro- 
dynamic significance. Recalling that 2x[y¥(z, p) — (zx, 0)] represents the total flow 
through the circle of radius p centered on the z-axis and lying in a plane perpendicular 
to the axis, we have within an unessential additive constant 


vie, )= 2] f Suduao (4.1) 
where 
r= (2? + yp’ + B® — Qyub cos 6)” (4.2) 


denotes the distance of the point B(0, 6) in the plane @ = O from the point Q(z, u, @). 
The surface integral in (4.1) is the solid angle 2 subtended at B(O, b) by the foregoing 
circle whose trace in the meridional half-plane is P(x, p). We write, 





W(x, p) = 5- 20, b; x, p) (4.3) 
2r 

It is, however, more expedient to determine y(z, p) from (2.5) and (3.6). Accordingly, 

dy _ _ a6 _ m [” bp cos 6 — p* 
“= Op ot Jo R* a6 (4.4) 

so that 
m f[* [* bp cos 0 — p° 

Ay = v2, 0) -— 0,0) ==] | OS? dodr (4.5) 


Performing the elementary integration with respect to z, 


_m [* (bp cos 6 — p )x 
— = T / (b° + p’ — 2bp cos @)R aid (4.6) 








and using the transformations (3.3), (3.7), 





: »K _ »K an* ¥ 
—2mzp | pews 4mbzp(b p) sn° (u, k) du (4.7) 


oY = So + OF, Je mi(b +p) Jo (b+ p) —kri sn’ (u, k) 


With a view toward transforming the last integral in (4.7) into the standard form of 
the elliptic integral of the third kind, we introduce the complex parameter a through 
p—b 


sn (a, k) = ss en (a, k) = x. dn(a, k) = > (4.8) 
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and reach 


one . ; 
am rp ‘ i] . 
Ay =- : K + —TI(K, a) 
T (p + b)r, 2 


ok 2 4 
k>sna ena dnasn’ u 
; du 


where 


(Ix, a) = 


1 — k°sn°asn’ wu 


is the standard form of the complete elliptic integral of the third kind. 


(4.9) 


(4.10) 


From (4.8), sn (a, k) > 1 which characterizes the integral as being of circular type 
according to Legendre’s classification. Here an explicitly real representation in terms of 


a real parameter ¢ is obtained by putting 
a= K+ ‘te 


— 2K <e< 2K’ 


The integral (4.10) can now be expressed as follows: 


(K, a) = I(K, K + Ze) 
= KZ(K + we, k) + nz 
= KE(K + ze, k) — (K + iE + nm 
valid within the open ranges 
(2n — 1)K’ < € < (2n + 1)K’ 
in which » denotes an integer. For the values 


e = (2n + 1)K’ we have n(K. a) = 0 


(4.11)° 


(4.12) 


(4.13) 


(4.14) 


Z(a, k) and E(a, k) stand for Jacobi’s Zeta function and the incomplete elliptic integral 
of the second kind of the argument a. By addition theorems and Jacobi’s imaginary 


transformation for sn u, en u, dn u, E(u) as well as Legendre’s relation 
KE’ + EK’ — KK’ = 
equations (4.8) to (4.13) lead to 


¥(2, p) = — E ee a ple K | 
T Tr 


with 
A(e, k’) = KE(e, k’) + (E — Kye = KZ(e, k’) + = 
2K 
and 
sn (e, k’) = —, en (e, k’) = ie dn(e, k’) = et? 
Ts 2 Ps 


(4.15) 


(4.16) 


(4.17) 


(4.18) 


6K’ and E’ designate the complete elliptic integrals of the first and second kind respectively for the 


modulus k’. 
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The appearance of n in (4.12) and its disappearance from (4.16) and (4.17) is ex- 
plained as follows: The integral (4.10), having as path of integration the straight line 
segment (0, K), is single-valued in the neighborhood of x = 0, p = b. The quantity a 
has the form given in (4.11), a = K + de, where e increases by 4K’ in a circuit around 
(0, b), and II(K, K + ie) is periodic in ¢ with the period 2K’. On the other hand, 
KZ(K + ze, k) is infinitely many-valued around (0, b). Now, in the relations (4.12), n 
has to appear to produce single-valued right-hand members, the left-hand member being 
single-valued. However, since y varies continuously in a circuit around (0, b), it is 
correct to replace II(K, a) in the process of computation of y by KZ(a, k) plus a single 
additive constant for the entire range —2K’ < « < 2K’. That is what has been done 
to obtain the results as expressed by (4.16) and (4.17) from which n has disappeared.’ 

5. Discussion. Comparing (4.18) and (3.17), we note that 


y = am(e, k’) 

5.1 
a [ - dt = (5 ) 
~ Jo (1 — (k’)’ sin’ t)'” 
where F(y, k’) is the incomplete elliptic integral of the first kind. The functions k’(2, p) 
defined by (3.13), (3.5) and e(x, p) given by (4.18) or equivalently by (3.17) and (5.1) 
may be regarded as curvilinear coordinates in the half-plane p > 0. The coordinate 


lines k’ constant (k = constant) are the “iso-modular’’ circles 
vr)? Je _ 1+k)? ? _ 4k’)? . 
(7) + E k* 7 k* 6% 


\s any iso-modular circle is traversed in a counterclockwise sense, ¢ increases by 4K’. 
In view of (5.1), 


e=0 for y¥=0 


e= +K’ for v +5 


f (5.3) 
e= +2K’ for y=cuAn 





lime = y 
k’—0 
so that « tends toward the polar angle y as P(x, p) — B(0, b) (see Fig. 1). 
The auxiliary function A(e, k’), defined by (4.15), is of basic significance for what 
follows. A(e, k’) is multiple valued for —o < ¢ <o. For a closed path encircling B 
once we obtain for the circulation of A 


A(e + 4K’, k’) — A(e, k’) = 2 (5.4) 


Here we face the alternative of either admitting a many-valued A and permitting a 
circulation, or of cutting the semi-plane along x = 0 0 < p < Db so that A becomes 
single-valued but discontinuous along the cut. The ranges for e and y as previously 
announced by (4.11) and (3.17) correspond to such a cut. 


7The arbitrary additive constant was chosen consistent with ¥(0, p) = —m for p > b. 
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The local behavior of A in the vicinity of point B is characterized by 


lim A(e, k’) = y (5.5) 


k’—0 


Confining the discussion to the range —2K’ < e« < 2K’, we record the values: 





r>0,p=0 1=9(2 -1) 

“ro 
2<0,p=0 1=9(2 +1) 

ary 

(5.6) 

z= +0, p < 6 A=-7 
r= -0O,p<b A=f 
S20. 9 > 6 A=0O0 


Equations (5.6) reveal the discontinuity of A along « = 0,0 < p < b in presence of a 
cut. 

It follows from (5.5) and (5.6) that in absence of a cut the stream function y(z, p), 
given by (4.16), is also a multiple valued function whose circulation is twice the total 


strength m of the source ring. 





W(e + 4K‘ k’) — Ye, k’) = 2m (5.7) 
Again, for the cut domain, —2K’ < e« < 2K’, we find in further agreement with Wein- 
stein® 
(—2m for xz>0 
y(z, 0) =4 
| 0 for a | 
¥(+0, p) = —2m for p<b (5.8) 
¥(—0, p) = 0 for p<b 
| 
¥(0, p) = —m for p>ob ) 
Observing 
lim y = m(? _ 1) (5.9) 
T 


f—a—0 


we recognize the local significance of y in the neighborhood of P as essentially that of 
the polar angle y. 

As pointed out earlier, the multiple-valued or, alternatively, discontinuous character 
of ¥(z, p) was to be anticipated in the presence of a multiply connected domain of regu- 
larity of ¢(z, p) whose boundary here consists of the z-axis and the point B(0, b). It 
should be pointed out that the cyclic nature of the stream function is also immediately 
evident from the solid angle interpretation underlying (4.3). 





8See loc. cit. 
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Figure 2 shows the streamline pattern of the source ring, i.e., the curves ¥(z, p) = 
constant, for the quadrant x < 0, p > 0 of the meridional half-plane. The curves given 
correspond to equal increments in y. The diagram also shows the locus v,(x, p) = 0.' 

6. Results for the source disc. We consider next a homogeneous distribution of 
sources of total strength m over the circular region x = 0,0 < p < b. The determination 
of the velocity potential ¢(z, p) and of the associated stream function ¥(z, p) may be 
achieved by integration with respect to b of the corresponding functions for the source 
ring. We cite directly the results obtained. 
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arb” 


—2m 
= —— 
(x, 0) lzl+fe 
¢- == as rz — 0 
rb 
—2 : 
v.(2, p) = — oe E K + 4 
rb r; 
. _ __+2m ; a 
v(x, 0) = ro(To x lz |) for M 0 
2m : 
v.(+0, p) = +7 for p<b 
) 
v(0, p) = 0 for p>b 
it 
v.—- ae as rz - 0 
rb 
v,(z, p) = am | e+ 3 ua K — | 
ab’ p ri 


v(x, 0) = 0 


v, — om | to 7. 2 — 0 as r2 > 0 
rb is 


m | , x(x? + 2p? + 2b’) 
p) = - ak - = 


rb ? 


V(x, 


(—2m for xr>0 
W(x, 0) =) 





0 for t <0 


¥(+0, p) = m( 2, _ 2) for p<b 


¥(—0, p) = -—S for psb 


¥(0, p) = —m for p>b 


y—> —™m 





—K+()° — p yA — rh’ 
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(6.3) 


(6.4) 
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In agreement with Weinstein, we note that Stoke’s stream function is single valued 
in the simply connected domain of definition whose boundary consists of the entire 
z-axis and of the cut along the segment x = 0, 0 < p < b. Moreover, ¢, v, , v, , and y 
are single-valued throughout the foregoing region which corresponds to 0 < k’ < 1, 
—2K’ < « < 2K’. We observe that v, is constant over each one of the two faces of the 


disc and that 


n(+0, p) — v(-0, 9) =F =o, O< psd (6.5) 


where o is the flux per unit area of the source disc. The velocity component v, , on the 
other hand, has no discontinuity as we pass from one face of the disc to the other, but 
has a logarithmic singularity at point B, i.e., at the edge of the disc. These results are 
in agreement with the general theory of source sheets. 

The solution corresponding to a homogeneous source distribution over an annular 
region bounded by two concentric circles is of course at once obtainable from the results 
just given by application of the principle of superposition. 

7. The vortex ring. We now turn to the flow generated by a vortex filament of cir- 
culation T along the circle x = 0, p = b. According to a theorem of vortex theory,” the 
velocity potential at any point of the space due to a single closed vortex filament is 
proportional to the solid angle which it subtends at that point. Thus for the case under 
consideration, using the notation introduced in (4.3), 


d(x, p) = == Q(x, p; 0, b) (7.1) 
4r 


Comparing (4.3) and (7.1), we note a reciprocal relationship between the stream function 
of the source ring and the potential of the vortex ring. By virtue of (7.1), (4.3) and 
(4.16), the explicit formula for ¢(z, p) is immediate. 

Alternatively, we recall'® that the flow corresponding to a closed vortex filament is 
identical with that induced by a uniform distribution of doublets over any surface 
bounded by it. The axes of the doublets are assumed to be everywhere perpendicular to 
the surface, and the surface density (strength per unit area) of the doublet sheet is equal 
to I'/4z, where T is the circulation of the vortex. The solution for the vortex ring there- 
fore coincides with that for a homogeneous doublet disc, which in turn is obtainable by 
differentiating with respect to x the solution for the uniform source disc given in the 


preceding section. Thus, 


‘ ‘ —bT a a s 
¢ (Vortex Ring) = ri aq? Source disc) (7.1) 
ete. 


Carrying out the required differentiations, we reach: 


‘See, for example, H. Lamb, Hydrodynamics, 5th ed., University Press, Cambridge, 1924, p. 195. 
See H. Lamb, loc. cit., p. 195. 
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o(2, p) = 2 E K + N 
2r Lr; 


o(x, 0) = #2 (1 lz ) as 2a2>0 or 2«<0 








To 
’ wr A 9 
¢(+0, p) = Fs for p<b (7.2) 
¢(0, p) = 0 for p>b 
ly 
—> as > 0 
on d 7 J 
rh. . v-2#- 
v(x, p) =>5 jx +e p B| 
2nr; ? 
| 
Y 2 
v(x, 0) = +5 t (7.3) 
“ro 
| 
r | cosy d 8r, a ) 
v, + | a (10g 7, sin 7 — 0 as r, > 0 
Tz iz F b° . 
v(2, p) = z . E 2 2 =< K - K | 
“7 pr; To | 
v,(z, 0) = v,(0, p) = 0 f (7.4) 
lr sin y 
v.— - ie. — 0 as To > 0 
riz 4 i 
Wx, p) =5 2 Tg 7B 
2r ; 
; 
¥(z, 0) = 0 (7.5) 
y Sr 
y— a2 to oD is 2| — 0 as r, > 0 | 
2r Te 


The formula for ¥(z, p) in (7.5) coincides with that given by Maxwell."' In contrast 
to the results for the source ring, the stream function here remains single valued if the 
cut is removed whereas the potential ¢(z, p) in (7.2) becomes cyclic having a circulation 
equal to the circulation of the vortex ring: 

o(e + 4K’, k’) — oe, k’) = T (7.6) 
Equation (7.6) follows also directly from (7.1) and agrees with the general theory of 
vortex filaments. 


11S8ee H. Lamb, loc. cit., p. 219. The cortesponding streamline pattern appears on p. 221. 
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8. The vortex disc. We consider, finally, a distribution of concentric circular vortex 
filaments over the circular region z = 0,0 < p < bin which the radial circulation density 


be given by 


daw ae (8.1) 


h? 


so that the circulation around the annular region bounded by the circumference of the 
disc and the circle of radius p is 


C(p) = | c(t) dt = rf — (2) | (8.2) 
: ; b 
The total circulation of this disc-vortex is 


CO) =P (8.3) 


II 


It follows from a theorem on vortex sheets’® that the foregoing vortex disc is equivalent 
to a non-uniform distribution of doublets over the same circular region. The axes of the 
doublets are again parallel to the z-axis and the surface density (strength per unit area) 
of the doublet sheet is C(p)/4r-. 

The solution for the vortex disc under consideration, however, is obtained most con- 
veniently by integration with respect to b of the solution for the vortex ring discussed 
in the preceding section. We merely state the results of these lengthy computations. 


' 
g(x, p) =; 5] 3ar,E — 4 (x? + 4p’)K + (227 — p’? + D*)A 
2rb rT; 
4 ! 1\2 
d(x, 0) = F 5p (ro — | x |) as > 0 or z<g 
, (8.4) 
I'r,72 
¢(+0, p) = —s for p<b 
2b 
d(x, p) = O for p>b 
¢—0 as rT, > 0 ) 
r 20” ) 
v(t, p) = - .| ane —-"- K+ 2ra | 
rb rT; 
7 — 
v(x, 0) = Br (ro — | x |) r (8.5) 
, 8 , 
= (4 — log ai — 0 as ro > 0 





2See, for example, Handbuch der Physik, vol. 7, Julius Springer, Berlin, 1927, p. 43. 
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lr" | ar, S43 ) <= 
v(t, p) = 3 K — — (x b° 2p" — pl 
A(X, p aE _ + b° + 2p°)K oa | 
v(x, 0) = 0 | 
; r : 
v(+0,p)= +55 for p<b * (8.6) 
) 
; 
v,(0, p) = 0 for p>b | 
v, + ry > 0 as r, 2 O 
, rb 
I TiN" eo K 20 —2x — db’ 
W(x, p) 3 ( Ea +bhp — »') - aa c : rik + 2p | | 
rh 3 1 3 | 
| 
Yr, 0) = 0 > (8.7) 
2Tl | 
yo as r, > 0 
or ) 


It is seen that ¢, v, , v, , and y are single-valued in the simply connected region 
0 < k’ < 1, —2K’ < e < 2K’. In contrast to the solution for the source disc, v, has a 
logarithmic singularity at the edge of the disc but is continuous along x = 0,0 < p < b. 
The velocity component v, has a finite jump discontinuity as the disc is traversed. By 
(8.6) and (8.2), and in agreement with the theory of vortex sheets, we obtain 


2M%p_  __ dp) 


v,(+0, p) — v,(—0, p) = =3 O0<p<b (8.8) 
b dp 


9. Remark on numerical evaluations of results. The velocity potentials, velocity 
components, and stream functions occurring in the four solutions given in this paper, 
involve beyond elementary functions exclusively the complete elliptic integrals K and 
{ and the cyclic function A(e, k’). In view of (4.15), (4.16) and (5.1), the values of A 
are readily obtained by aid of tables of complete and incomplete elliptic integrals of 
the first and second kind. Alternatively, tables of Jacobi’s Zeta function and of Jacobian 
elliptic functions may be used without recourse to tabulations of incomplete elliptic 
integrals. An extensive numberical tabulation of A values has been completed in the 
process of work now in progress at Illinois Institute of Technology. 
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STUDIES ON TWO-DIMENSIONAL TRANSONIC FLOWS OF A 
COMPRESSIBLE FLUID.—PART II* 


BY 
S. TOMOTIKA anp K. TAMADA 
University of Kyoto, Japan 


7. An alternative method of treatment. In the foregoing sections, it has been shown 
that the fundamental non-linear equation in the ¢, y plane governing the flow of our 
hypothetical gas has exact solutions of practical interest in several cases, and for two 
of these cases, the corresponding flow patterns have been discussed in detail. Unfortu- 
nately, however, it seems very difficult to proceed further along similar lines of attack, 
since it is unlikely that any more exact solutions of practical interest can be found. 

In Part II, therefore, an alternative method of treatment is developed on the basis 
of the fundamental equation in the hodograph plane. Just as in the case of the isentropic 
flow of the real gas, the fundamental differential equation governing the flow of our 
hypothetical gas becomes linear in the hodograph plane, and by suitable linear com- 
binations of appropriate exact solutions of this linear differential equation, various flow 
patterns of practical interest have been obtained. 

If the independent and dependent variables be interchanged in Eqs. (3.5), we obtain 
the fundamental equations of motion for our hypothetical gas in the w, @ plane in the 


forms: 

pu = 2kwys ’ De ” Vw ’ (7.1) 
where, as before, 

eo _ ¥+i - 

w = | od, kat. (7.2) 


Eliminating yy from these equations, we get a partial differential equation for de- 
termining the velocity potential @ of the form: 


Pu ae . Pu _ 2hwoss = 0, (7.3) 


while elimination of @ yields an equation for determining the stream function y of the 
form: 
V,, oo 2hwwoo = 0. (7.4) 
The flow of our hypothetical gas can be found by solving either of the above two 
equations. Since, however, the latter equation is more simple in form than the former, 
discussions on the flow of our hypothetical gas will be made on the basis of the latter. 
We shall use in place of w a new variable z, defined as 
z= v/2kw. (7. 
Equation (7.4) assumes the form 


Vi. — Zoo = O. (7.6) 


d Jan. 21, 1949. Part I of this paper appeared in this Quarterly 7, 381-397 (1950). 


* Lece 1vé 
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This is a simple typical partial differential equation of the “mixed type’’, changing from 
the elliptic to the hyperbolic type according as z < 0 (i.e. w < 0) or z > 0 (i.e. w > 0); 
in other words, according as the flow is subsonic or supersonic. 

The equations of the characteristic curves for equation (7.6) are easily found to be 

0— 0 = +3”, (7.7) 

where @ is an arbitrary parameter. 

8. Two different families of particular solutions of Eq. (7.6). To solve our funda- 
mental equation (7.6), we first assume y in the form: 

y = Zze"'”’, (8.1) 
where » is a constant. Then, putting this in (7.6), we get an ordinary differential equation 
for determining the function Z(z) of the form: 

dZ = 

i +y2zZ = 0. (8.2) 

dz 


The solutions of this equation can be expressed in terms of Bessel functions of orders 
+ 1/3 and thus we obtain the following family of particular solutions of equation (7.6), 
namely: 

Z rads (202""*) | ¢ “ 


y=: x 4 (8.3) 


} ivd 


On the other hand, another family of particular solutions of equation (7.6) can be 
obtained in the following manner. 
We now put’ 


z =, 6 = nz," (8.4) 
and transform the independent variables in (7.6) from z, 6 to & ». Then, equation (7.6) 
becomes 


9 » 15 ‘ 
Eve: — 3tnv:, + | > 1) Yon + r ny, = 0. (8.9) 


To solve this equation, if we assume y in the form: 
y = 2 Y(n), (8.6) 


with an arbitrary constant r, we get the following ordinary differential equation for 
determining the function Y(n), namely: 


9 d°Y 15 F dY : 7 
(2a - 1) 7+ — 3r)n +r(r—1)Y = 0, (8.7) 


1This transformation is suggested by Eqs. (7.7) for the characteristic curves. It may be remarked that 
the curves defined by » = +2/3 coincide with the two characteristics passing through the point (z = 0, 
6 = 0). 
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and this equation is further reduced to the form: 


ad. ae i ~{t-3 \epa¥ 2 — 
yt Be det lo 3 3ar}s dk g M(r 1)Y = 0, (8.8) 


when use is made of a substitution of the form: 


9 « ; 
{7 = 6. (8.9) 
Equation (8.8) is a special case of the hypergeometric equation and therefore its 
solutions can be expressed in terms of hypergeometric functions. Thus, we find that in 
case when | ¢| > 1, 


- qi ya-r) pf 1 jum 5_1,,4;1) 
| § r(4 "), & girs ea 


Y = 2 (8.10) 
(1/3)r _ l 1 saul I 2 . 1) 
| § r 3 r, 2 33 ’ c , 


while, in case when | 1 — ¢| < 1, 


(8.11) 


| 
slike TS Oh ee a ) 
—_ ‘- ! d = s)r iY) = o = ne = . — 
j|a- 9) fe re ‘): 
Hence, we ultimately obtain another family of particular solutions of our funda- 
mental equation (7.6) in the following forms: for | @| > (2/3) 2°”, 
5 1 4. 4) 


g (2/3) (1l—r) ( ee = > one. 
= F\3 (lL — 1), § — 3°53 5 og 
) 
| 


y= . (8.12) 
oon 505-375 Soe) 
and for |6@| < [2(2)'*/3]2"”, 
[er(ha - 9, - $8 - Frit - 9), 
= 4 cesieouth - : (8.13)? 
(MY bade ded Bet = 29) 


Furthermore, these solutions can be extended beyond the respective limits by the use 
of the well-known formula for analytic continuation of the hypergeometric function. 
In this place, the relation between the two families of particular solutions (8.3) and 


2In the range (2/3)z2 < |@| < [2(2)"/*/3]z*/*, the second expression becomes complex. But, it is 


readily found that both its real and imaginary parts are still particular solutions of equation (7.6). 








130 8S. TOMOTIKA AND K. TAMADA [Vol. VIII, No. 2 


(8.12) will be considered. For this purpose, we remember a formula’ for a definite integral 
containing Bessel function of the form: 


ax se ny ; bh? 
| e° "J \( bv)" dy «( 3 Pa (a ’ Atett ; A +> i; *). 


2a a“ ( 


“0 


If, in this formula, we take 


- 1 a oe 

A= $3, #=37 3% 
_ 2 3/2 | 

a= 8, b=3 | 


and multiply both sides by z'”, we readily find that the right-hand side of the formula 
becomes, except for a certain constant, just equal to the expression on the right-hand 
side of (8.12). 

Therefore, the expressions for y as given by (8.12) can also be put in the following 


forms: 


paz’ . , duo wr rn ae (8.14) 


so far as the integrals on the right-hand side are convergent. 

Comparing this expression with (8.3), we can easily understand in what manner the 
expressions for y as given by (8.12) are constructed by superposition of the former family 
of particular solutions (8.3). 

9. Properties of singular points of the solutions. The solutions defined in the fore- 
going sections have been found independently by other investigators [1], [2]. However, 
we may now generalize these functions in such a way as to render them more useful. We 
note that the fundamental equation (7.6) retains its original form when the variable 6 
is replaced by +(@ — 7d), where J is a real constant, and therefore, the expressions as 
obtained from (8.12) and (8.13) by replacing @ by +(@ — 2A) are also solutions of equa- 
tion (7.6). As will be seen presently, this transformation of variables, though of a very 
simple nature, yields very important results. 

If this transformation of variables is applied to the second expression in (8.13), we 


obtain the expression for y in the form: 


( a/ ~ \2) (1/6) + (2/3)r oy ‘ — 
J 9(8 — 2d)" \ l l | l ( 2 9(@ — 2A) 
y = 241 — > r( +-—-y, be r-l— (9.1 
v 12° J aes * Ss" 673 4z ) 
Since this expression contains the factor: 
9(6 — 7 2) (1/6) + (2/3)r 
{1 - =) ; (9.2) 
42 ) 


the function y will have a point of singularity at the place where this factor vanishes, 
the order of singularity being determined by the constant r. Since, however, the value 


3G. N. Watson, A treatise on the theory of Bessel functions, Cambridge Univ. Press, 1922, p. 385. 
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of r can be chosen arbitrarily, we can make the above expression (9.1) for y to have a 
point of singularity of any order at the place where the relation: 
42° — 9(6 — ir)” = 0 (9.3) 


holds. 
Remembering that z and @ are both real, the above relation (9.3) can be split up into 


the following two relations: 


yA 4 .¢ 
2+ : (A — 6) = 0, 1 3 NO = 0. (9.4) 
When A = 0, these relations give 
9 ; 1/3 
6=0, z= -(2 »’) , (9.5) 


Thus, in this case there is an isolated point of singularity on the axis @ = 0 in the range 
z < 0, i.e. in the subsonic region.* 

On the other hand, when A = 0 the points of singularity of the function y as given 
by (9.1) are distributed continuously along the curves: 


which are nothing else than the two characteristic curves passing through the point 
(0 = 0,z = 0). 

Thus, summarizing the above results we may infer that the isolated point of singu- 
larity of any order in the subsonic region in the hodograph plane undergoes remarkable 
changes as it approaches the supersonic region, and as we pass from the subsonic to the 
supersonic region, this concentrated singularity is prolonged to the continuous distribu- 
tion of singularities along the above two particular characteristic curves in the supersonic 
region. 

[t will naturally be expected that such an interesting behaviour of singularities would 
still be found in the case of solutions of an equation of the mixed type of the following 
form: 

Vir + (Cz + C22" + +++) hoo = 0, 
where ¢; , @ , -** are all constants, which is nothing but an extension of equation (7.6). 

10. Special cases in which the solutions are expressible in terms of elementary 
functions. In this section, we shall consider a few special cases in which the hyper- 
geometric functions involved in (8.12) degenerate into elementary functions. 

First, we take, in (8.12), 


(_ 3 a 
| 2’ 3, 
r= or = 
1 
— 53 —2. 


‘The position of the point of singularity is definitely determined by the value of \. Physically speak- 
ing, the constant \ corresponds to a parameter for determining the Mach number of the uniform flow 


in the problem we shall consider. 
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Then, replacing 6 by +(@ — 7d) as before and remembering the well-known relation: 
F(n, 8; 8; xz) = (1 — 2)”, 


we obtain the results that 











12(@ — id) ae , : 2 ; x... a) = < orn 
6’°3’3’9(@— ay) f uw i ‘se 
y= 4 (10.1) 
9 ~ 
(6 — id) F(A ic i= — =- — 
6’3’3 ’°9(6 — id)’ | beg ae \ 
4$(@ — 7A) — =e 
+) 
(00 — ay F [2s yon ) a! 
S° 8 7S ° OO — aX)" f ae ae haan 
(6 — aa)? — =a 
8) 
y = (10.2) 
; er) ae ae 4z 6 — ir 
(A — > (: a eS eames ’ = — = — 
| i) F\3 6 +3 * 90 — a? 


3 . ae i‘. 7/6 
[ he = wy - 4a} 


It is worth noticing here that the two expressions in (10.2) can evidently be derived 
by differentiating respectively the two expressions in (10.1) with respect to 6. Similar 
results hold good in general; namely, as is suggested by the form of the fundamental 
equation (7.6)-as well as by the form of its general solution as given by (9.1), we can 
obtain, by differentiating or integrating any solution with respect to @, another solution 
whose order of singularity differs by unity from that of the original solution. 

In the next place, we take 

r= 1/2 


in (8.12). Then, replacing 6 by +(@ — 7A) and remembering the general formula: 


F(n, n+552n+ 1; 2) Z (2) (4/1 — x —1}”, 
2 Me 
we have 


( , 
ew awl £22445 ..08 
‘oie Ait + gt 
| 


g\'/3 0 mE we 1/3 
= (2) f (@— 1A)" — 4 +(0- a} 








Y 

y= 4 (10.3) 
jen ave. t okay t,o 8 aaa) 
| (6 — a) F( 6? 672? 67 1396 — =a)’ 


» 
~ 
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and by differentiation or integration of these solutions with respect to 0, we can also 
derive, as before, various solutions in closed forms having different orders of singularity. 
It will readily be seen that all the above solutions (10.1), (10.2) and (10.3) have, in the 
case when \ =~ 0, an isolated point of singularity at the point @ = 0, z = —(9/4 *)” 
in the subsonic region, while in the case when \ = 0, their singularities are distributed 
continuously along the two characteristic curves 6 = +2/3 2°”. 

11. Solutions having logarithmic singularity. When the parameter r takes the value 
of —1/4, the two solutions in (8.13) coincide with each other. In this case, there exists 
another independent solution which has logarithmic singularity and represents the flow 
due to a source or a vortex in the hodograph plane. 

To obtain such a solution, we shall start again from the differential equation (8.8). 
When | 1 — ¢| < 1, the two independent solutions of this equation are given by (8.11), 
but in case when r = —1/4 these two solutions coincide with each other. Another inde- 
pendent solution can be obtained in the usual way. Thus, making use of the fundamental 
solutions of the form (8.11), namely: 


= oe a Oe eae ) 
y.= FAG r), 36 3h l c), 


, (1/6) + (2/3) r ry 1 1 1 1 7 2 
ili MAR t gn gt ag t grt — tp 


which are valid in case when r # —1/4, the required solution Y, can be obtained as: 
y 3 . Y, _. 7 
oo 35 me Te le’ 
where the factor 3/2 has been multiplied in order to regulate the form of the solution. 
Calculating the limiting value, we ultimately obtain, in case when r = —1/4, the 
following independent solutions valid when | 1 — ¢| < 1: 


ris = F . a Tee = - _ 
V,(s) Me aya :) 1+ Dial — 9", 
1 2 


Y.(¢) = Yi log 1 -HO+ La Dv + T mes - 9 (1D 
n=1 m=0 m + DP m+ 75 


1 5/5 a «tte 1 a ) 
a = 3 (5+ 1) (; “ 1) bi + y (Gtr A} 
Hence, in case when r = —1/4, the stream function y for the flow due to a source in 


the hodograph plane can be obtained by the use of the second solution Y,(f) in (11.1), 
taking (8.4), (8.6) and (8.9) into account. We thus have 











pe eS oe 9(6 — ay") ( 9(6 — ny’) 
tities AS. gti - si 
— cf 43 bf ( _ 96 — ay} 
+ 2a 2 5 t [-aeire ae > (11.2) 
+ig ™tT ip 


where a, assumes the same value as that given in (11.1). 
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Further, by differentiating this expression with respect to @ once or several times, 
we can obtain, as before, the stream function for the flows due to a doublet and a multi- 
plet having logarithmic singularity of any order in the hodograph plane. 

[t is worth mentioning here that logarithmic singularities of these solutions experience 
also the characteristic change as mentioned in §9 when we pass from the subsonic to the 
supersonic region. 

12. An example of flow. We shall now construct a field of flow by employing various 


solutions obtained in the foregoing sections. 
In the first place, by superposing the first two solutions in (10.3), we have 


, Bk mt 
v {/(o — Ca) 9 w —(0- A) f 
j 4)? Sk oa oe 
— \a/ (4 “tay = 9 w+ (d— wed | , (12.1) 


where the variable z has been replaced again by the variable w with the aid of (7.5) for 


the sake of comparison with the analysis in Part I. 


6 








ee ee eee oe eee 
— ee ae 


Fic. 17. The flow represented by the stream function y . 


Taking the real part of the above expression (12.1) for y, we define a solution: 


vi, = Ry). (12.2) 


Then, this solution has evidently a branch-point of order 1/2 at the point @ = 0, 


w = —(9/8 k)'”°)’*” and becomes always zero along the axis 6 = O in the range 


w< —(9/8 k)'“*n*”. 
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Next, differentiating y in (12.1) partially with respect to 6, we have 


. 





( Kor _ )1/3 } = ; Bk 1/3 
11'S 6a; = 9 w —(0—a)f + \\ (@—tr)"— 9 w'+(0-2)} 
Wo = - = ——__—_——— (12.3) 
3 nso. Ons 
yaa) ma 3 w 
and, as mentioned already, this becomes also a particular solution of equation (7.6). 
Thus, by taking the imaginary part of this expression, we define a solution: 
v2 = I(W), (12.4) 
which has a branch-point of order —1/2 at the point 6 = 0, w = —(9/8 k)'”*r*”, just 


as in the case of the previous solution y, , and becomes always zero along the axis 6 = 0 
in the range w < —(9/8 k)'*”. 


19 
| 
| 








Fic. 18. The flow represented by the stream function y. 


Figs. 17 and 18 show schematically the flow patterns in the w, 6 plane as represented 
by the above two stream functions y, and y, respectively. 

These figures suggest that if we superpose these two solutions in the form ¥, + Ky, 
where K is a constant, we can obtain the flow pattern in the w, @ plane as shown in 
Fig. 19 (a), which corresponds, in the physical plane, to the field of flow past a sym- 
metrical body having sharp leading and trailing edges as shown in Fig. 19 (b). 

In such a flow, the condition of nearly uniform velocity is considered to be fulfilled, 
and therefore, the flow of our hypothetical gas is expected to approximate satisfactorily 
the corresponding flow of the real gas subject to the adiabatic law. 

Detailed investigations will make it clear that the flow can still be continuous and 
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irrotational in spite of the existence of limited supersonic regions in the vicinity of the 
surface of the body (Fig. 19 (b)), but at a certain critical Mach number there appear 
two points of singularity / = 0, one each on the upper and lower surfaces, which grow 
up into the curves of singularity J = 0 for still higher Mach number, thus violating 
the one-valuedness of the physical field of flow. However, we shall not enter into the 











a eo 
q<C i i. 
a % 
; €>C \ 


a —$§ 




















(b) Physical plane. 


(a) Hodograph plane. 


Fic. 19. 


numerical computations for the present example of flow, since it has been found that 
more full discussions concerning the similar field of flow can be made by introducing a 
second new hypothetical gas which represents the real gas to a higher degree of ap- 
proximation than our first hypothetical gas used in Parts I and II. The flow of our 
second hypothetical gas will be discussed in detail in Part III. 
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A NEW APPROACH TO PROBLEMS IN TWO DIMENSIONAL FLOW* 


BY 
MONROE H. MARTIN 
University of Maryland 


1. Introduction. The dynamical equations and the equation of continuity 
p(u,u + uy) + pz = O, plv.ut+ ov) +p, = 0, (pu). + (pv), a, (1) 


for the velocity components u, v, the density p and the pressure p in the steady two- 
dimensional flow of a fluid, subject to no external forces, constitute an underdetermined 
system of three partial differential equations for four unknown functions 


u = u(x, y), v = v(x, y), p = p(x, y) #9, p = p(x, y), (2) 


of the independent variables 2, y. 

To make the problem determinate, an additional relation may be prescribed, as is 
frequently done when it is assumed the density p is a function of the pressure p. In 
this paper we prefer, at the outset, to deal with the underdetermined problem (1), re- 
serving the right, however, to make the problem a determinate one whenever it is con- 
venient to do so. 

A novel feature in the investigation is the use of the stream function » and the 
pressure p as independent variables. With this device, the integration of (1), barring 
exceptional cases,’ is equivalent to the integration of a system 

htm =, kot wm =90, Gg = av, D), (3) 
of two partial differential equations for two unknown functions — = &(y, p), 7 = n(y, p). 
Corresponding to the indeterminacy in (1), the fluid velocity g = q(y, p) in (3) is an 
arbitrary function” whose specification characterizes the system (3). 
Once a solution of (3) has been found, the flow is presented parametrically by 


i —n ly, P); = Ev(y, P); gé=¢- pt, ’ 
(4) 
y = &(y, p), v=m¥,P), 2=17-Pm- 
For y const. the equations in the first and second columns constitute parametric 


equations respectively for the stream lines in the physical plane and their hodographs 
in the hodograph plane, the pressure p serving for parameter. In the third column &, 7 
denote lift and drag functions. 

Provided the flow does not have a degenerate hodograph, the variables y, p may be 
eliminated from (4), and the usual transformation from the (u, v)-plane (hodograph 
plane) to the (2, y)-plane (physical plane) is obtained. 


*Received March 24, 1949. Presented September 9, 1948 at the Madison meeting of the American 
Mathematical Society. This investigation was carried out under project NOL 139(Tm5-1) “Rotational 
Gas Flows” of the Naval Ordnance Laboratory, White Oak, Silver Springs 19, Maryland. 

1See Theorem 1 in Sec. 2. 

2The phrase “arbitrary function” is used in a restricted sense throughout the paper to designate a 
function subject to such continuity restrictions as may be necessary to insure the validity of any argument 
in which it enters, but which is otherwise entirely arbitrary. 
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The form of equation (3) is invariant under the change of parameter 
y= vy"), 


a property that leads immediately to the Munk-Prim substitution principle’ for the 
special case of two dimensional flows. 

By eliminating one of the unknown functions in (3), a Monge-Ampére partial differ- 
ential equation is obtained for the other. In view of the symmetry of (3), the unknowns 
£, n satisfy the same Monge-Ampére equation. 

The application of the Cauchy-Kowalewsky existence and uniqueness theorem to 
this Monge-Ampére equation for a prescribed regular analytic function q°(y, p) leads 
immediately to an existence and uniqueness theorem for flows having a preassigned 
analytic pressure and velocity distribution along a prescribed analytic are C. The func- 
tions £, 7 in (4) are presented by power series 


t=h+hi(W—wW) +-::, n=ko +h’ — Wo) +--:-, (5) 


in which the h, k& are regular analytic functions of p which are readily calculated from 
the initial data by elementary algebraic operations and differentiations. 

A necessary and sufficient condition for (3) to have a solution (5) linear in y — po 
is that the flow be irrotational. For a solution of this type the flow defined by (4) is the 
well known Prandtl-Meyer flow,* and a simple geometric construction results for the 


stream lines. 

The study of flows for which é, 9 are polynomials in Y — yp offers an attractive 
problem for investigation left untouched in the present paper and should lead to in- 
teresting examples of rotational gas flows. 

If the drag function é (lift function 7) defined in (4) be introduced in place of £(n) 
as unknown function, and y, y(y, x) taken for independent variables by an Ampére 
contact transformation, the Monge-Ampére equation for é(n) is replaced by a quasi- 
linear equation for (7). This leads to a general formulation of the method followed by 


Kiebel” in his investigation of rotational flows. Unfortunately, a simple parametric 


representation for the flow, analogous to (4), no longer appears to be available. 

In the final paragraph the treatment of one-dimensional unsteady flows is reduced 
to a minor formal modification of the method followed in the study of two-dimensional 
steady flows. 

2. The equations of motion. As is well known, the system (1) can be written in the 


alternative form 


(Dp + pu’), + (pu), 0 (puv), + (p ~ pu’), = (). (pu), + (pv), = 0, 


3M. M. Munk and R. C. Prim, On the multiplicity of steady gas flows having the same streamline pattern. 
Proc. Nat. Ac. Sci., 33, 137-141 (1947). The principle is established by Munk and Prim for three dimen- 
sional flows. 

‘T. Meyer, Ueber zweidimensionale Bewegungsvorgaenge in einem Gas, das mit Ueberschallgeschwin- 
troemt, Dissertation, Goettingen (1908) or Forschungsheft 62, V.D.I., Berlin (1908), pp. 31-67. 


digkeit s 
A. Kiebel, An example of an exact solution of a two-dimensional rotational flow in gas dynamics 


) & 
(Russian). Prikl. Mat. Mekh. 11, 193-198 (1947). English translation by E. Rabkin, Technical Trans- 
lation No. TT-35, National Research Council of Canada, Division of Mechanical Engineering. Ottawa 


(1947). 
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from which it follows’ that , given a solution (2) of (1), there exist three functions £, 7, 
of x, y such that u, v, p, p, &, 7, ¥ constitute a solution to the system of Pfaff 


dt = puv dx — (p + pu’) dy, dn = (p + po’) dx — pu dy, dy = pu dx — pudy. (6) 


The variable y is, of course, the stream function, and to lend the variables £, 7 phy- 
sical significance, we write (6) in the form 


dé = udy — pdy, dn =vdy + pdxz, dy = pudzx — pudy. (7) 
Along a streamline 
di = —pdy, dn = pdx, dy = 0, 


from which the physical meaning of £, 7 is apparent. The difference between the values 
of £ at two points on a streamline forming the boundary of a solid immersed in the fluid 
equals the horizontal component of the force experienced by the body due to fluid 
pressure along the intervening are. Since 7 may be interpreted as the vertical component 
of this force, the functions &, 7 are frequently referred® to as the drag and lift functions, 


respectively. 
If we introduce auxiliary functions 
f=E+py, n=1-— pe, (8) 
of x, y, the system of Pfaff (6) becomes 


dé cdy + y dp, dn = vdy — x dp, dy = pv dx — pu dy, (9) 


. | 


which suggests that ¥, p be taken as independent variables. This restricts us to those 
solutions (2) of (1) for which no functional relation exists between the functions y, p of 
x, y. In other words we rule out at this point those flows for which p = p(y), that is, 
flows for which the streamlines are isobars. Such flows form a very restricted class which 
will be studied in Sec. 3. 

With the introduction of ¥, p as independent variables in place of x, y, the functions 
E mn, p of x, y become functions of y, p 


& = &(y, p), n= ny, p), p = ply, p) # 9, (10) 


with x, y, u, v, —, 7 given as functions of y, p by (4). Substituting for x, y, u, v from (4) 
in the last equation of (9) and remembering that y, p are independent variables, we 
find the following conditions on the three functions (10) 


Evtyp t+ mmr te =90, yt + nym, = 0. (11) 


Corresponding to a solution (2) of the underdetermined system (1), we are led by the 
above considerations to a solution (10) of the underdetermined system (11). 

Let us consider the converse problem. Given a solution (10) of the underdetermined 
system (11) can a solution (2) of the underdetermined system (1) be constructed? 


‘See, for example, Bateman, H., The lift and drag functions for an elastic fluid in two-dimensional 
irrotational flow. Proc. Nat. Ac. Sci., 24, 246-251 (1938). The notation in this paper differs from ours, 


with Bateman writing —y¥, —X, —Y for what we have designated by y, &, 7. 
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The underdetermined system (11) may be regarded as a determinate system of two 
partial differential equations for two unknown functions &, y, the function 
p = ply, p) = O, (12) 
being regarded as an arbitrary function whose specification characterizes the system. 
From this standpoint, the converse problem is formulated as follows. The function 
(12) being specified, given a solution £, 7 of the determinate system (11), can a solution 
(2) of the underdetermined system (1) be constructed? 
If we assume that 
o(x, y) 
O(y, p) 
the equations in the first column of (4) determine y, p as functions of x, y. Substitution 


for , p in the equations of the second column, and in (12), then yields u, v, p as functions 
of xz, y. We shall now show that the four functions (2) of x, y so obtained constitute a 


= Ey aMvp bi Ny vévr = J # 0, (13) 


I 


solution of (1). 


For this we need the following formulae derivable from (4): 
Juz = Eyutoo — Ein » Tut, = E59M — Every » 
Jv, = ae re a Sun%es ? Jv, = Nyy lpp — Nip , 
(14) 
Jp , EppPy wae EypPp ’ J p, = NopPY —~ NerPd » 
Jp. = —typ; Jp, = —Ny» 
and the formula 
Py = P(Eyn + Myo — Eveton — 290%), (15) 


obtained by differentiating (11) partially with respect to p, y. On substituting in (1) 
from (14), it follows from (4), (11), and (15) that (2) is a solution of (1). 

When (13) and the second equation in (11) are solved simultaneously for £,, , np» , 
it is found that 
(16) 


Ep = pd ny , Np» = —pJé, , 


Pp 
and consequently, from (4), that 
pqJ =2,+y,, q =u +v, 
so that the condition 


z, + ¥ > 0, 


and (13) are equivalent as long as p ¥ 0, gq ¥ 0. The last condition amounts to requiring 
that the directional derivative of the pressure along a streamline remains finite. 

The above results may be summed up in the theorem. 

TuEorEM 1. The solution of the underdetermined system (1) and the solution of the 
determinate system (11), containing the arbitrary function p = p(y, p) are equivalent 
proble ms provide d the streamlines are not isobars, the directional derivative of the pressure 


along the streamlines remains finite, the density does not vanish, and the flow has no stag- 


nation points. 
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The specification of the density function p = p(y, p) amounts to assuming that along 
each streamline the density is a function of the pressure alone, the functional dependence 
being permitted to vary from one streamline to another. 

Returning to the first equation in (11), we may use (4) to write it in the form uu, + 
vv, +p = 0. Holding y fast, we integrate this equation from pp to p to obtain Bernoulle’s 
theorem 

l 1. a 

50 =50—] e'dp, m= PlW), d= aol¥)- (17) 

2 2 Fn 
Here po is some conveniently chosen reference pressure, not necessarily the same for all 
streamlines, and q, is the associated fluid velocity, both po(W), go(¥) being arbitrary 
functions. 

Thus the system (11) is replaced by the equivalent system (3) in which the fluid 
velocity g = qg(¥, p), given by Bernoulli’s theorem, characterizes the system (3). Apart 
from the condition gq, ~ 0, the function qg may be taken arbitrarily, with 


p = (qq), do = Wy, polW)).- (18) 
On differentiating (17) twice with respect to p, it will be found that 
M* — 1 = qq Qo ; (19) 


where 1 = q¢/G is the Mach number. The flow is accordingly supersonic or subsonic as 
a O or i... < ws. 
For po = const and a density function of the form p = f(p)/¥(y) we find 


l a? l @ 7 v(y)P(p) Go = q (y) P(p) ™ [ & ° (17’) 
2° 2 Jno S(P) 
In the special case of a perfect gas f(p) = p", Viv) = k" with 0 < n < 1, we put pp = 0 
and find 
2 l 2 e l—n witht 
9 q = 9 = er ) . (17 ) 
Here 3q5 is the total head, k = k(W) = exp (S/c,), with S denoting the entropy and c, 
the specific heat at constant volume. For an incompressible fluid with variable density 
f(p) 1 and 


Ll; 
Gq =5% — V(Y)(p — po) (17’”’) 


-_ 


is a linear function of p. 
The formulae (14) lead to an interesting expression for the vorticity wo = u, — v,. 
Direct substitution from them, followed by the elimination of £,, , n,, by (16), leads to 


= = —sle (20) 
For a perfect gas with uniform total head this reduces to 
n adSp 


o= gee oe 


l1—ndye,’ 
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to imply a theorem due to Crocco,’ which states that the vorticity is proportional to 
the pressure along a streamline. 

The expressions for the lift and drag functions, given in (4), follow at once from (8), 
in view of the first column in (4). 

3. The exceptional case p = p(w). As has been pointed out above, the flows for which 
the streamlines are isobars must be treated separately. According to Prim,” in the case 
of a perfect gas, the streamlines for such flows are either concentric circles or parallel 
straight lines. We shall show how a modification of our previous analysis generalizes 
Prim’s result to any fluid for which the density is a function of the pressure along a 
streamline. 

Under this assumption, p = p(y) implies p = p(y). Let us take y, s for independent 
variables, where s denotes the are length measured along a streamline from some fixed 
point upon it. One finds that (9) is replaced by 


dit = (u+ yp’) dy, dn = (v— xp’) dy, dl = pvdx — pudy, ’=d/dy. (21) 


/ 


Thus £&, » are functions of y alone and we may W rite 


, ; 


- ; ; 1 “oer 
£&=ut yp, n =vU— ap, vty — Uy =p , vz, — uy, = 0, (22) 


the last two equations coming from the last equation in (21), since y, s are independent 
variables. From the first two equations we find 


+ yp’ = 0, v, — x,p’ = 0, (23) 


u 
8 | 


which, provided p’ + 0, may be used to eliminate x, , y, from the last equation in (22) 


to obtain wu, + vv, = 0, and consequently 


u" + y= q; q = ql). (24) 


When uw, v are eliminated from (24) with the aid of (22), we obtain 


r\2 eh\2 
- ( 
(: + 1) — fi — =; ) = 4 (25) 
p \ p 
This shows that any streamline is a circle with center (—’/p’, &’/p’) and radius 
r=q p’ =r y ). (26) 
To prove that the circular streamlines are concentric, we employ the identity 
(Ey (2) P 
7 “ir 7 y= a — 2 . 
Pp \p 
Since the right side is a function of y alone and wu, v depend on s, this implies 
¢/ n’ 
=; = const., 7 = const., (27) 
P p 


7Crocco, L., Eine neue Stromfunktion fuer die Erforschung der Bewegung der Gase mit Rotation. 
Zeits. angew. Math. Mech., 17, 4 (1937). 

®R. C. Prim, On the equations of plane rotational flou of a perfect gas in natural coordinates and on plane 
rotational qas flow with constant velocity magnitude or constant vorticity along the streamlines. Naval Ord- 


nance Laboratory Memorandum 9264, (1947) SB Project No. 19.2. 
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and, as a corollary, that 


Il 


qr’ —p' =0. 
If we write (26), (27) in the form 


I 
oP» 


y’ = @, qr’ = 1 
dr q py F 


and eliminate 7’, we obtain 
dp 2 
"dr 9 
from which it is clear that any two of p, g, p may be specified functions of 7, and the 
third is then uniquely determined, if we except p which is only determined up to an 
additive constant. 
If p’ = 0, it follows from (22) that the streamlines form a family of parallel straight 
lines. Each streamline is a straight line along which u = ~’, v = 7’ are constant. To see 
that these lines are parallel, let 


r= 2% + 8 cos 8, Y = Yo + ssin 8, (28) 


where x, , ¥% , @ depend only on y, be the parametric equations of the family of stream- 
lines. Substitution in the third equation of (22) yields 
q(x,sin 6 — y) cos 0) — g0’s = p'," 


which, since p p(w), implies 


9’ = 0, pq = (xj sin 6 — y/, cos 6)", (29) 
so that 6 = @, = const., as stated, and the Cartesian equation of the family of stream- 
lines is 

x sin 0 — y cos & = n, (30) 


where » denotes the length of the normal let fall from the origin to a streamline. If we 
take 2) , Yo in (28) to be the coordinates of the foot of this normal, we have 


n’ = xjsin 0 — ys cos %, 
from (30), so that the second equation in (29) reduces to 
pqn’ = 1. 


The following theorem has accordingly been proved. 

THEOREM 2. For fluids in which the density is a function of the pressure along a stream- 
line, of the streamlines are isobars, the streamlines are either concentric circles (p # const.), 
or are parallel straight lines (p = const.). In the first case, any two of p, q, p may be pre- 
scribed as functions of the radius r to determine the third as a function of r. In the second 
case, p, g may be taken as arbitrary functions of the distance n of the straight line from the 
origin. 

4. The Munk-Prim substitution principle. A change of parameter Yy = ¥(y¥*) along 
the isobars leaves the form of the system (3) unchanged, the new system being 


we 2 *2 _ #2 7 a 
Sve + 7 =F » as ae + V*Npp = 0, 
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where 


( b(y*), DP); n* = n(y(y*), P); q* = = Qu(y*), P)- 


dy - 

ee ee ee q ee : + 

w= —y = —, ee ie Soe, OCU UE pe = §, 
ay 

tees ek ¢ «ea oe — oY ~ ae ‘+. = 

P= C= & FH = Tews "9% Pap = 2, 
¥ 


so that the new flow has the same streamlines, pressure distribution, lift and drag 
functions, as the old system. On the other hand, we have 


dy dy ( $+) 4 
re a i cs 
’ 2 dy* ”? p que) ” 


a u 
dy* 
the last relation being a consequence of Bernoulli’s Theorem (17). 
Thus we are led to the following substitution principle first formulated by Munk 


u 


and Prim” for steady three-dimensional flows. 

From a given flow, a new flow having the same streamlines and pressure distribution 
may be derived by multiplying the magnitude s of the velocity vectors tangent to a given stream- 
line by the same factor \ and the density by the factor \”~, the factor d being permitted to 
change from streamline to streamline. 

By properly choosing the factor d either go(W) or Y(y) in 
constant. Thus, as Munk and Prim have pointed out, in the case of a perfect gas the 
new flow may be taken to have either constant total head or constant entropy throughout. 

5. The Monge-Ampére partial differential equation. The problem of integrating (3) is 
equivalent to the solution of a Monge-Ampére partial differential equation. To see this, 


obtaining 


17’) may be reduced to a 


we solve (3) simultaneously for ny , 7), 


and consequently 


9Un — SySy 
; p2 ; 
Ss 


2 i72 Iv et gpawten 
\q = ) \q - 24) 
from which a Monge-Ampére partial differential equation for & = &(y, p) 
20d vEyn — Wvéveor TF Iv — (99 + Wey + 7 Evveoo — Evo 


arises when 7 is eliminated by partial differentiation. 
Given a solution ~ of (33), we determine 7, up to an arbitrary constant by the 


S 


line integral 


= | iNny AW + 1» Ap}, 34) 





9Loc. cit. 
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with the aid of (32), and then determine 7 by the line integral 


7 = | iny dW + n, dp} 


with the help of (31), (34). From the symmetry of (3) it is clear that 7 is also a solution 
of (33). 

We propose to investigate the existence and uniqueness of flows when the fluid 
velocity g in (33) is a prescribed function of y, p. From (18) this amounts to prescribing 
the associated fluid velocity qo = qo(W) and the density function p = p(y, p). 


Consider an are 


e: x= xp), y = y(p), >» SPS His 
in the physical plane. Can this arc be taken to be one of the streamlines of a flow in which 
the pressure of the fluid at the point (x(p), y(p)) equals the parameter p? Furthermore, if 
such a flow exists, 1s it uniquely determined? 
For C to be an are of a streamline as required, it is necessary and sufficient, from (4), 


that a constant y exist such that 


(Wo, P) = y(p), Ey(Wo , P) = (bo, p) cos O(p), (35) 


itr 


where 
6 = are tan y’/2’, ’ = d/dp, 


is the angle of elevation of the tangent line to C. If (35) is replaced by 


tivo, p) =et+ | y(p) dp, Ey(yo , P) = Q(vo, p) cos Op), c = const., (36) 


these conditions give the initial data for Cauchy’s problem for the partial differential 
equation (33). To apply the Cauchy-Kowalewski existence theorem,'” we write (33) in 
the form 


Ey, = [Ey — 2k vev> + WvEeeo — VP Iov + (Ie + IEG Ev) (37) 


[t follows from this theorem that if g(¥, p) is regular analytic for | y — yo! < 6, pm S 
p < p, and if C is an analytie are with x” + y” > 0 holding for p, S p S p, , there 
exists exactly one analytic solution £ of (37) meeting the initial conditions (36), provided 
that 

qv) , p) # 0, two, Pp) = y #0 for Py Spsp.- (38) 
This solution permits the expansion 


=h + hi(W — Wo) + °° thy — wo)" +-:-, lv —wwo| <b < 6, (39) 


¢ 
Ss 


valid for sufficiently small 6, in which h, = h,(p) are regular analytic functions of p for 
MW =PEP 

The function n(y, p) is now determined by (34) after substituting from (39) for & in 
(32), the arbitrary constant in (34) being adjusted so that ,(Wo , po) = —2(po)- 


Hadamard, J., Lectures on Cauchy's problem in linear partial differential equations, New Haven 


(1923) pp. 9-16, in particular p. 16. 
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The restriction y’ ~ 0 in (38) is not an essential one. If per chance y’ = 0 at some 
point of C, the réles of £, 7 are interchanged in the analysis and (38) replaced by 


avo, p) € 0, Np(Wo, p) = —2’ € O. (40) 


By hypothesis x’, y’ do not vanish simultaneously on C, so if (38) fails, we can always 
be assured that (40) holds. 

The condition g(¥> , p) ¥ 0 is equivalent to requiring that there be no stagnation 
points on C. 

The above results are summarized in the theorem below. 

TueoreM 3. If C: x = 2(p), y = y(p) is an analytic are with x” + y” > 0 holding 


for py. < p S p., of the fluid velocity q(y, p) is regular analytic for | ~ — Yo| < 4, Po < 


p Sp, with qo , p) 0, the arc C is a streamline for exactly one flow (of preassigned 
sense on C) with pressure p and fluid velocity q(Wo , p) at the point (x(p), y(p)). The stream- 
lines and the hodographs of the flow in which C is imbedded are presented parametrically 
by (4) in which the functions é(p, p), n(y, p) permit the expansions (5) valid for sufficiently 
small | ~ — |, the coefficients h, k being regular analytic functions of p in the interval 
Do Sp Sp- 

To calculate the coefficients in (5), these expansions are substituted in (3). On the 


assumption 


g = q +a(y—y Ss qo ~ Y, $1) 
where q@ , 9: **: are regular analytic functions of p for pp» S p S p, , a routine computa- 


tation shows that the h, k are determined by 


> (r + 1)(n + 1 — r)(h,. Nass + k w 3 -r) = Ga» 


42) 
YH a ea + ee has) = 0. 
In particular, for n = 0, 1, 2, we find 
w+ =q, hi'h, + kik, = 0, 
, 1 ‘aA Ad l Ad ‘7 
hike + kyky =, hhothg + ki’ k2 = — 5 (hi'hy + Ki'h,), 
4 2 
13) 
Jo 2 
Lie + Eb ow Be SG + I 
6 3 
aa - _ - sade “ 
hé’hs + ko’ks = — 3 (ho’h, + ki’k,) -— 3 (hi h. + ky’'k.). 


The coefficients h, , ky are, from (4), determined by 
hs = y(p), ki = —2x(p). 


The first pair of equations in (43) are solved simultaneously to determine h, , k, and 
the second pair may then be solved simultaneously for h, , k, , since the determinant 


isk; — Rit,” = |¢ x’ + y’?)]! . 
ig y )} 


1950] A NEW APPROACH TO PROBLEMS IN TWO DIMENSIONAL FLOW 147 


does not vanish by hypothesis. Proceeding in this manner, the remaining coefficients 
are determined step by step. 

6. Prandtl-Meyer flows. Leaving aside the case where the density or velocity 
vanishes, it is clear from (20) that a necessary and sufficient condition for the flow to be 


irrotational is g, = 0. 
For (3) to have a solution (5) 
—E= ho +h (wv — yo), n=k+hkh Ww — Wo), (44) 


linear in ¥ — yp, it is clearly necessary that g, = 0, that is, that the flow be irrotational. 
Conversely, we shall now show that among the irrotational flows there are solutions 
(44) of (3). 


We set 
q = 4(p) (45) 
in (3) and take q,., = h, = k, = 0 for n = 2, 3, --+ in (42), or (43), to get following set 
hii+k =q, hyht’ + k, ki’ = 0, hho’ + k,ko’ = 0, (46) 


of necessary and sufficient conditions for (44) to be a solution of (3). To satisfy the first 


pair of equations, we place 


h, = q cos 8, k, = qsin 8, (47) 
and find that 
Dp rr\ 1/2 
6 = Op) = %& + [ (") dp, (48) 


to determine h, , k, in (46) as functions of p. 
Leaving aside, for the moment, the determination of hy , ko , we substitute in (4) 
from (44) to obtain 


x= —ki — 8'(W — Wo) cos (0 — py), u 


q(p) cos 6(p) = u(p), 


(49) 


y= hi —s(¥ — yY)sin (6 — yp), v q(p) sin 6(p) v(p), 


where g(p), @(p) are the functions of p mentioned in (45), (48), and s denotes the arc 
length of the hodograph curve 


H: u = u(p), v = v(p), 
in the hodograph plane. The angle u defined by 
cot un = g6’/q’ (50) 


is the angle between the radius vector OP to a point P of H and the perpendicular ON 
let fall on the tangent line to H at P (see figure below). From (48) and (19) we find 


cot » = (M?— 1), 8! = 8'(p) = —(Gp)", (51) 


so that » = u(p) is identified as the Mach angle. 

On placing p = const. in (49), it is clear that the isobars are straight lines along 
which all velocity vectors have the same magnitude and make the same angle » with 
the isobar. These straight lines are the “‘Mach lines” in the flow defined by (49). 
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With ¥ = const. in (49), the streamlines are presented parametrically in terms of 


the pressure p by the equations in the first column. From the second column it is clear 


that all streamlines have the same hodograph, a characteristic feature of Prandtl-Meyer 


flows. 


We return to the determination of hy , ky . Taking Y = Wp in (49), we obtain an initial 


streamline 
x= 2(p) = -—k, y=yp=h, 


in terms of which the last condition in (46) becomes 


y’ cos 6 — x’ sin 6 = 0, 52) 
on substituting from (47). Thus we may take any curve 
C: x= f(r), y = g(7), (at) + (da) > 0, 
dr dr 
with continuous non vanishing curvature for the initial streamline, for if we write 
up) = f(r), Y\p) = gr), T = TP) 
and introduce the slope m = m(r) of C, condition (52) becomes 
m(r) = tan 6(p) 53) 
an implicit equation for the determination of the parameterization t = 7(p). In other 


words, once C is given, the pressure distribution along it is determined by (53). 


The above considerations are illuminated by the following geometrical construction. 




















Fig. 1. 


At a point 7 with parameter 7 on C in the physical plane we draw the tangent line 
TU. From O draw OP in the hodograph plane parallel to 7U to meet the hodograph 
curve H at P with parameter p. This determines geometrically the correspondence 
7 = 1(p) between the parameter 7 at 7’ and the parameter p at P affected analytically 


by (53). 
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The isobar through 7 is constructed by drawing a line 7'J parallel to the normal ON 
let fall upon the tangent line to H at P. To construct a streamline, the length (Gp) is 
laid off from 7 along 7J to determine a point S. The locus of S as T traces out C is a 
streamline. The remaining streamlines are “homothetic transformations” of this stream- 
line, obtained by laying off the lengths (¥ — ¥o)(pG)~* along TI from T. 

7. The drag and lift functions. The Monge-Ampére equation (33) is carried over into 


a quasi-linear equation 
Evy + 2qGE Ev» — (0° Aer — (940 + WEE — Ibe = 9, 9 = av, -&), (54) 
with y, y in place of ¥, p as independent variables, and the drag function & = ty, y) 
as unknown function, by the Ampére contact transformation 

v=¥, D=y, E=E-y, t=u F=-p. 
With y¥, y serving as independent variables, (7) is written 


dé =udy — pdy, dy = prx, dy + p(vx, — u) dy, 


S 


from which we see that 





u=f,, p= —t,, t, = (gv), zt, =w'. 
Since u° + v° = aq’, the last two equations may, from (18), be written 
. 99» . gy 
vy eC 1/2? ty 7 42\1/2? (55) 
ig — &) ig — & 


from which x may be eliminated to obtain a second derivation of (54). 
Here p is presumed to be replaced by —é, . Given a solution £ = &(y, y) of (54), the 


Sy 


Cartesian equation of the streamlines is 
r=2x(y,y) = | {zy dy + x, dy}, 


where x, , 2, are the functions of y, y given in (55). The pressure p and velocity com- 


ponents are, of course, given by 


Formulas analogous to the above have been obtained by Kiebel with —y, x as inde- 
pendent variables for the special case of a perfect gas, without, however, explicitly 
computing (54) or noticing the physical significance of ¢. To obtain the actual formulas 
of Kiebel, the lift function 7% is used in place of the drag function é, the lift function 7 
becoming the function’ x used by Kiebel in his investigation of flows for which x has 
the special form 

x= —H(yp)e- 1-4, 
8. Unsteady one-dimensional flow. If we disregard the second equation in (1) and set 


v=l yh 


NOp. cit. p. 2. 
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in the remaining equations, we obtain the differential equations 
plu,u 4- U,) + P: = 0. (pu), aa Pp = 0, 


for unsteady one-dimensional flow. We find (9) replaced by 


dé udp + t dp, dy = pdx — pudl, 56) 
and taking y, p for independent variables, * 
we find, instead of (11) 

Egy, —ty tp = 9, £,e,, — % = YU, 57) 
where p = p(y, p) is an arbitrary function. Eliminating x by partial differentiation, the 
role of (33) is taken over by 

l “i 
Evuspp — Syn = , 8) 
\P/ p 


a differential equation which was obtained by another method by the author’ in the 


special case of a perfect gas. 
Given a solution § = &y, p) of (58), the flow is presented by 


z=2(¥, p) = / (Ey 


in view of (56), (57); or, provided £,, * 0, so the last equation can be solved for p, the 


+ p') dy + &&,, dp}, u = E,(Y, p), t= §,(¥, DP), 


Jr 


equations defining the flow can be put in the form 

z= zy, i), u= uly, t), p = ply, t). 
For fixed y these equations give the positions x, velocity u, and pressure p for a gas 
particle as functions of the time t. 


2The case p = p(y) is disregarded. 


Martin, M. H., The rectilinear motion of a gas II, Amer. J. Math. 57, 410 (1945). 
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WALL EFFECTS IN CAVITY FLOW—I* 


BY 
G. BIRKHOFF, M. PLESSET, anp N. SIMMONS 
Harvard University; Naval Ordnance Test Station, Pasadena; Ministry of Supply 


1. Introduction. We consider below high-speed cavity flow of a liquid past a solid, in 
which a long cavity, filled with gas (air or water vapor), is formed behind the solid. For 
such a flow, the cavity pressure p, is presumably nearly constant, and the cavitation 


parameter Q or K is defined asf 
Q = K = (p,; — p.)/4 pv; - (1) 


Here p, is the free stream pressure, and v, the free stream velocity. There is considerable 
experimental evidence’ that for small K (say K < 0.15), the physical conditions assumed 
in classical “wake” theory are approximately fulfilled. We present various theoretical 
considerations, which indicate strongly that a free jet is preferable to an ordinary water 
tunnel with fixed walls in studying such flows, because of two related wall effects. 

First, we deduce in Sec. 2 the existence of a blocking constant for each tunnel and 
model, below which A cannot fall. 

Second, we show that a large wall correction must be made, for drag coefficients 
made in a water tunnel with fixed walls. The case of “infinitely long’’ cavities is covered 
in Secs. 7-8 below; the case of “finite cavities” will be dealt with in Part IT. 

The numerical results presented in Secs. 7-8 can be deduced from formulas of Réthy. 
In fact, some similar results were obtained by Valcovici [5] in 1913; but our results are 
more exact, much more extensive, and based on simpler computations. They suggest 
methods of estimating wall corrections, entirely analogous to those suggested by Prandtl 
and Valeovici [5] for ordinary flows with wake. However, we have felt it necessary to 
give a fresh discussion (in Sec. 9), which will be physically reliable for cavity flows with 
small K (say K < 0.2). It is notorious that the wake interpretation is entirely inaccurate 
physically. 

The formulas of Réthy have been extensively generalized by Mises [3]. We give, in 
Secs. 3-5, a further generalization, which permits one to determine any flow with free 
streamlines, whose hodograph is a circular sector. In Sec. 6, we discuss a new method for 
basing effective numerical computations on these formulas. 

2. Blocking constant. Consider the idealized cavity flow in a water tunnel with fixed 
walls, depicted in Fig. 1. We suppose an incompressible, non-viscous liquid, and a sta- 
tionary liquid-gas interface, with negligible turbulence. Further, we suppose a uniform 
upstream flow with velocity v) , and a uniform free downstream velocity v, as the cavity 
approaches its maximum cross-section A, , in a tunnel of cross-section Ay . The rate 


*Received April 8, 1949. The material of Part I was developed by Birkhoff and Plesset in 1947 (see 
Abstract 54-7-258t, Bull. Amer. Math. Soc., 1948). The material of Part II was developed by Simmons 
about the same time (see Proc. 7th Internat. Congr. App]. Mech., London, 1948, vol. 2, p. 601), and the 
concept of a “blocking constant”’ was introduced by him. 

tIn Part II, we shall use Q to avoid confusion with elliptic function notation. 

‘Mostly unpublished, cf. P. Eisenberg and H. L. Pond, Water tunnel investigations of steady state 
cavities, David Taylor Model Basin Report 668 (1948). Also G. Birkhoff, Recent progress in free boundary 
theory, Proc. 7th Int. Congress Applied Mech., London, 1948, p. 7. 
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of increase of liquid momentum per unit time is clearly pvjA, — pvj>A4o , where A, = 
A, — A, . The total thrust on a long section of liquid is (p — p.)Ao — D where pp is 


the upstream pressure in the free stream, and D is the drag. Hence 


pv; A, — pto>Ao = (po — p.)Ao — D. 


But by Bernoulli’s equation, pp, + pvi/2 = p. + pv;/2; by conservation of volume, 
v,A, = A> . Substituting, 
I 2 2 l 2 9 
D = pAgds (Xi — %) + Wo — uo) = 5 pA,(v; — vp). (2) 
By Bernoulli’s equation, applied to (1) with py = po, 
K = (v,/y)° — 1, or (1+ K)'” =20,/y. (2’) 


Hence if we define the drag coefficient, as is usual, in terms of the upstream velocity, so 
that 


Cp = 2D/mA = (:/v9 — 1)°Ad/A, (3) 

we have the exact relation 
Cp = (14+ K)”’ — 1)A,/A S K°A,/4A. (4) 
This defines a blocking constant, analogous to that occurring in sonic flow with J/ near 

unity: 

Keay "s/Ay”. (4’) 
Since C, ordinarily varies between .0625 and 1.00, we have in practice (A,/A)'”*/2 S 
Kymin & 2(A,/A)'”’. Thus, to achieve K = .05, we must have A,/A about 400, at least. 


3. Mathematical assumptions. Henceforth we shall discuss theoretical plane flows 
with free streamlines, usually (cf. [2], Chap. XI) called flows with “wakes’’. For greater 





Q 








physical realism, we shall refer to cavities instead of wakes. In Part I, we shall consider 
mainly infinite cavities behind wedges in the center of tunnels with fized walls (Fig. 1), 
free jets (Fig. 1a), and bounded jets issuing from orifices (Fig. 6). 

In addition to the drag coefficient defined by (3), for such flows, we shall consider 
the drag coefficient 


C, = 2D/p?A = Cy/(1 + K) (3’) 
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based on the downstream velocity. Following Valcovici, we show that C, gives a wall 
correction which is much smaller than C, . In fact, we show that for wedges, the correc- 
tion is infinitely smaller. 

Since the drag coefficient is in theory independent of dimensions for a given value 
of a/b, one need only consider the case in which the downstream velocity is unity; 








further one may put a = z/v where v denotes the upstream velocity. Thus, the stream 
function is normalized so that it goes from zero on BOSB’ to x on CC’. In this notation, 


since a = m/v, 
x (1 —v) = 
Co == ae eine 9 
b v (5) 
——— T (1 — v)’ -? 
C,= b (5’) 


Similarly, one finds from consideration of momentum and continuity for the open 


tunnel, or free jet, case 


a, 
C = 2-(1 — cosa) 
b 
and, if the previous convention regarding the stream function is adopted, a = /v, 
v = 1 so that 
Y T a) 
C= 25 (1 — cosa). (6) 
) 


It is to be noted that the upper half of the flow may be regarded as a jet issuing 
from an angular orifice. In the fixed wall case, fhe ratio v of upstream velocity to down- 
stream velocity is simply the “coefficient of contraction”. Numerical correlations be- 
tween v and b for this case, and between a and 6b for the case of a free jet, have been 
obtained by Mises [3] who obtained excellent agreement with experimental observations 
on jets. 

4. Technique of conformal transformation. Analytical formulas covering the cases for 
which we have derived numerical results, in Secs. 6-7, may be found in many places.” 


2M. Réthy, Klausenburger Berichte (1879) and [4]; U. Cisotti, Rendic. Palermo 28 (1909), 307-52 
and Idromeccanica piana, Milan, Secs. 140, 146; V. Valcovici [5]; R. Von Mises, [3] and [1]; P. Frank and 
R. von Mises, Partielle Differentialgleichungen der mathematische Physik, Ch. XI, Sec. 2. 
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However, we believe that our formulas provide a simpler basis for effective numerical 
computation than any in the literature; certainly, our numerical results are new. 
Let the physical plane be the z-plane with z = x + zy, and the origin at the stagnation 


point 0. The complex potential is W = U + iV, where U is the velocity potential and 

V is the stream function. The simplest mathematical representation is obtained by 

choosing conventions as to sign so that ¢ = dW/dz is the negative complex conjugate 
Rf - plane 


) B,C ciB' 








0 B,C 


£ + in of the vector velocity —£& + zy. In Figs. 2-2a we show the hodograph or ¢-plane 
for the fixed wall and free jet cases. The hodograph is a sector of the unit circle, of m/n 
radians. 

We shall now present a unified treatment, applicable to all cases in which the pre- 
ceding hypothesis is satisfied (and also if 1 ¢ is a circular sector). 

Whenever this is the case, the one-to-one conformal or “‘schlicht”’ transformations 
coe, ¢ 0-01 4+ 2)/01 — &), and ¢ > y, where 

y= (1+ °)°/1- ¢) (7) 

map the hodograph onto a half-circle, a quadrant, and the upper half-plane, respectively. 
The most general schlicht transformation doing this is therefore’ given by 


Avy+B A1+20° + "l ae (8) 


Cy + D Cl+ Mme + Ie" g 


Here A, B, C, D, i, uw are real, and X ¥ uy, since otherwise we would have ¢ — const. 


3C. Caratheodory, Conformal mapping, Cambridge, 1932. 
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Now in the case or Figs. l-la, the W-graph, or region of the W-plane corresponding 
to fluid points, is clearly the infinite strip 0 < V < . Hence e” occupies the upper 
half-plane; this is the case treated by Réthy and Mises (op. cit. supra). 

We now show that this restriction on the W-graph is unnecessary. The W-graph 
can be a plane, half-plane, or infinite strip, with or without cuts. The reason is simply that 
in these cases, the upper half r-plane can be mapped conformally onto the W-graph by 
a Schwarz-Christoffel* transformation of the special form 


dW = R(r) dr, (R(r) a rational function.) (9) 


rm . ” . . . . 5 
Thus in the case of a symmetric wedge in an infinite stream treated by Bobyleff,” the 
W-graph of half of the flow is the upper half-plane. In any such case, we have 


dW = dr = (qdp — pdq)/¢. (9.1) 


In all the cases of jets (many of which can also be interpreted as flows with cavity or 
‘“‘wake’’) treated by Réthy and Mises,” the W-graph is an infinite strip. Here we can set 


dW = dr/r = (qdp — pdq)/pq = dp/p — dq/q. (9.2) 


The W-graph is a cut plane in various cases, including flows past oblique plates (Ray- 
leigh) and asymmetric wedges.’ In any such case, we can normalize to 


dW = 7dr = p(qdp — pdq)/¢. (9.3) 


The W-graph is a cut half-plane in other interesting cases. These include simulations 


tev" 
T - plane 





T =o T=! T=s? 


SRS: Sea 
B' Ss e) B,C C' 
Fic. 3. 








of a seaplane float’ by a plate or wedge under a free surface. In such cases, we can 
normalize to 
dW = rdr/(r — 1) = p(qdp — pdgq)/¢q'(p — g) (9.4) 


This will be simplified further in Sec. 4, formula (12’). 

‘See for example [2], Ch. X. 

5Jour. Russ. phys.-chem. Ges. XIII (1881). 

*Refs. [1], [3], [4]; also ref. [2], Sees. 11.51, 11.53, and Ex. 5 of Ch. XI. This case also covers the useful 
case of a bend in a pipe with straight walls, meeting in an angle on the outside, and in a curve on the 
inside, so as to get constant pressure along the curve. Such a design should minimize the tendency to 
cavitation and flow separation on the inside of the bends; a cavitation-free expansion can be constructed 
similarly. 

7See [2], Chap. XII, Sec. 12.50, and Exs. 3, 8; Rayleigh, Collected papers, vol. I, p. 287. 

8See A. E. Green, Proc. Camb. Phil. Soc. 31 (1935), 589-603; 32 (1936), 67-85 and 248-52; 34 (1938), 
167-84. In these papers special cases falling under (9.4) and (9.5) are treated. See also [2], Sec. 12.3. 
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The W-graph is a cut infinite strip (and the hodograph a circular sector) in various 
mr ° ° ee . ° — 9 ° 
other cases. These include an infinite free jet divided by a flat or wedge, a wedge in a 


T - plane 
s= cot nd/e2 


T =-s@ T =0 T=1 








>. oe ee 
C cB S 0 B 
Fia. 3a. 


> 


stream bounded by parallel walls, one of which is parallel to a side of the wedge (Fig. 3), 
and the jet produced by two streams joined in a wedge (Fig. 3a). In any such case, we 
can normalize to 


dw = r drt —— l 2 dr 4 a a 
(7 + alr —a — 1) ata’ T+a : na 


, 1 
qdp — pdq | a 1 a | 
ata’ q\p + aq) q(p — a ‘@) ; 


If the rate of flux in the two branches is equal, even though the flow is not symmetric, 
a = 1 and we have by elementary integration W = (1/2) Log (7° — 1). 

5. Algebraic simplification. The algebraic simplification in the complex domain given 
by Mises (refs. of footnote 2) for the case (9.2), is possible in the general case (9), when- 
as a linear combination of terms of 


ever R(r) = R(p/q) can be represented explicitly’ 
the form 


Rr) = I(r — a;,,)/M(r — ay). 


It is not always easy to determine R(r) from the physical data of the problem. But 
there always exist, by the theory of Schwarz-Christoffel transformations, real constants 


such that 


R(r) = >> - Bi + Bo + >> Bir — a). (10) 


Ss = a 


Here the a; represent points at infinity on tubes across which the stream function jumps 
x8; , while the a; correspond to stagnation points. 

Once these real constants have been determined (one may in practice have to carry 
one or two arbitrary constants through all calculations), each tr — a; = 7, = pi /q,; can 
always be expressed as in (8). Moreover, since 

%See footnote 8 and [2], Chap. XII, Exs. 2, 4, 6; the symmetrical case can be reduced to (9.2). We 
have sketched figures for the cases not in the literature. 

From the point of view of effective computation, the explicitness of the factorization is important. 
It can be very tedious to determine numerically all the roots of a general polynomial, and to verify that 


the determination is exact enough for subsequent computations. 
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d(Ln(p;/q,)) = dp;/p; — dqi/4 


dr,/T; 


dr = dp/q — pdq/¢ - (1) 


7, dt; = dp;/q; — p; dq;/q5 


dW equals ¢""' dé times a linear combination of polynomials in ¢", divided by denomi- 
nators of the form (1 + 2y,¢" + ¢°")"“’, where m(i) < 3. Now applying the general 
formula, dz = dW/¢ , we get after repeated synthetic division by factors q; , g, and q; . 

THEOREM 1. Let there be given any flow whose hodograph is a circular sector and 
whose W-graph is a strip, half-plane, or plane, with or without cuts. Then z(¢) can be 
expressed as a sum of integrals of the form 


[ ta; + bt") dt/(1 + 2n,¢" + 2)". (12) 


Here the m(z) are positive integers and the real coefficients a; , b; , v; can be explicitly 
evaluated in terms of the A, uw of (8), by rational operations. 
In the case (9.4), the reduction is especially simple since p — q = (A — uw)f" and 


p= (p-Qry. It yields 


dz = 2nt”* dt ‘+ w+ ‘< + (A — uw) (I + 3") 





A+u 
2 + ; 
q q (A — wt 
In any case, the quadratic functions 1 + 2»,¢" + ¢°" which occur in the denominator 
of (12) can be factored, in the real or complex domain. Thus 


2 , 
+ (A — 5} (12") 


(“+v ye" +v") if ly|>1 
1+ D+ pe = 4. if ly} =1 (13) 
e +o"ue -<™™) if ly| <1 


Here v = vy + (v — 1)” and a = (1/n) Cos’ v respectively, so that the factorization 
is explicit. We conclude, using partial fractions again. 

THEOREM 2. In the complex domain, we can represent 2(¢) as a sum of constant 
multiples of integrals of the form 


[ ve"? d¢/(¢” — B,)”, (m = positive integer). (14) 


In the case (9.2) treated by Réthy-Mises, we always have m = 1. A slightly weaker 
result holds in general. Since 


i. _ | _ _ xm — De" er her dt 
(" — By)" ("- 6)" ~  ("—B)™*’ 


we get through repeated integration by parts, the following result. 


THEOREM 3. Under the hypothesis of Theorem 1, z(¢) can be expressed as a sum of 
constant multiples of terms of the form ¢’/(¢" — 8,)”, and of integrals of the form 








[ e azn - 8). (14) 
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6. Effective numerical computation. In the general case, the effective computation of 
2(¢) involves various questions, which we shall only touch on. By introducing new 
variables —w; = ¢"/8,; , we see that the evaluation of (14’) can be effectively accomplished 
by constructing a table of the incomplete beta functions 


. 


t 
| w dw/(1 + w) (15) 
in the complex domain. After considerable study, we are impelled to the conclusion 
that this is a more effective way to compute 2(¢), in the general case when the hodograph 
is a circular sector than the method found in the literature, for a general n. We can 
certainly state 
THEOREM 4. Under the hypothesis of Theorem 1, 2(¢) can be explicitly expressed 


in terms of elementary functions and incomplete beta functions. 





- : : ‘ ; 1/h_k . 
If n = h/k is rational, then the substitution ¢ = @;’"c" reduces the complex integrals 
(14’) to the form f o” da/(co" — 1), where y, h are integers. Moreover 
a’ de . w de “Sky +1) y a , 
—h | ' [= =, _— — = 20 Log (¢ —w). (16) 
J —o ‘ —wo k=l 


Hence we conclude (cf. Mises [1], pp. 821-2, for the case (9.2)). 

THEOREM 5. If n = h/k is rational in Theorem 1, then we can express z(¢) explicitly 
in terms of complex logarithms and rational functions of ¢' 

However, unless k is small, the effective computation of a flow field by this method 
will be extremely small. Hence Theorem 4 however elegant, does not adequately de- 
scribe the real problems of obtaining numerical results of physical interest. Thus phy- 
sically, in the problems of Sec. 2, one is given n, and wishes to correlate the asymptotic 
velocity with the ratio b/a = (plate width)/(channel width); not A or yu. It is most 
convenient to correlate these as follows. 


In the fixed wall case, = 1 when W = —&,e” = 0;¢ =v when We" =o. Com- 
paring with e” = p/q, we get g = (1 — ¢")* and p = (¢" — n") (¢" — v "). This gives, 
after reducing to partial fractions, 

tad ans i 2 whe “ 
a= a. ap - = — nn (17) 
tm ¢ ey | eg 
° ; " a . — Ww 
In the free jet case, § = 1 when e” =o and ¢ = e’* = B is complex when e” = 0. 
This gives similarly 
Qnt"* dé nB"e"~ dé no" * dé _ 
f=. = = Soe aie eeerers (17a) 
— l— si er 


The application of the incomplete beta function to this “fixed wall’’ case is relatively 

simple. Thus if we write vf = e*'’"r, etc., in (17), we see that we could compute b(v), 

corresponding in the z-plane to ¢ = e*'’" from a table of the real incomplete beta function 
aie 


ig 2 n l —Il/n f ) 
&,(y) = | tr“ dr/i+ 7) = | w /" dw/(1 + w). (18) 


lt Jo 


We could also obtain the velocity (hence pressure) distribution along the wedge, the 


tunnel wall, and the dividing streamline from the same table. The case of the free jet 


is less simple, and will be discussed in Sec. 7. 








| 
; 
' 
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In this connection, we note that the real definite integral ,(1) can be expressed" in 


terms of the logarithmic gamma function y(n) = I’(n)/T(n). Specifically, 
( 1 ) (! )| 
#,(1) = — l=] = — — =z] |. 19 
2n E 2n 2n “\2 2n (19) 
Hence, if v is near one (case of large channel), we need only compute ®,(v) — ®,(1) 
and ®,(1/v) — #,(1), which are presumably easy, in terms of series (cf. Sec. 8). 


7. Numerical results for plate. In case the wedge is a plate, n = 2 and the preceding 
remarks are unnecessary. We can integrate (17) easily, getting 


9 
2 = —4Tanh™ ¢ + 2v Tanh™ ve + “ Tanh™ (¢/v) (20) 


b = (xr/v) — r + 2 — 1/v) Tan 2, (21) 
in the fixed wall case. The drag coefficient Cy and C, of formulas (5), (5’) may now be 
correlated with values of b/a, using v as a parameter. Numerical results are given in 
Table 1, and the results are shown graphically in Fig. 4. The “wall correction” for a 
closed water tunnel may thus be estimated. 

The numerical results demonstrate clearly that basing the drag coefficient on the down- 
stream velocity leads to a smaller correction. This is the result of Prandtl and Valcovici [5], 
but based now on much more extensive evidence. This conclusion will be discussed 


further in Sec. 8. 








In the free jet’* case, we get similarly from (17a) with n = 2, 
2a | 2dy (at __dg =) 
az = a Nii aie % 1 te 
l $ l + $ l Bs + Bs (20a) 
1 ( dé dt 
e\l— et 1+ rar 
b = r(1 — Cosa) — 2 Sina Log tan (= - a) (21a) 


This expression together with the drag formula (6) gives the values of Table 2. These 
values are also shown graphically in Fig. 4. 

In the ‘fixed wall” case of the closed tunnel, the pressure at the tunnel wall is a 
quantity which may be readily measured; thereby, a further comparison with theory is 
made possible. If pp is the static pressure and v is the stream velocity as a great distance 
upstream from the lamina, then the static pressure p, where the stream velocity is u, 


is determined by the relation 


PoP _ ue 
(p")/2°° vw l. (22) 


“See [2], p. 315; in [1], p. 110, a different notation is used. The function (zr) is tabulated in vol. I 
of The British Association Tables, and many other places.—Using this result, we can correlate Cp with 
the case of a symmetrical (Bobyleff) or assymetrical wedge in an infinite stream, in the case that the 
stagnation point it at the leading edge. Existing tables of the incomplete beta function do not cover the 
range nee ded for the present problem. 

“Tn the free jet case, if n = 2h is an even integer, we can also slightly simplify the computations since 


+ af* + 1 — af? + ) oo 1+ CB — @)i* + *. 
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Positions on the tunnel wall are z = x« + im/v and one gets at once from (20), with 
sus i 
1+ u 1 + vu 1 utov 
x= —2 log — + »y log —— + — log —_, 
—u 1 — vu v “u—v 


2 4.8 
= —4tanh™'u + 2 tanh vu + =tanh"'v-, 
v u 
where the origin of z is at the stagnation point 0. The pressure coefficient (22) has been 
evaluated as a function of x for three values of v; the results are shown in Fig. 5, and are 
tabulated in Tables 3a, 3b, 3c. 
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8. Case of bounded jet. A case of particular physical interest is that of a semi-infinite 
closed channel terminating in a free jet. The physical plane is depicted in Fig. 6, where 











FIXED C 
TUNNEL WALL 
~L_ 4 Vv 
| qd 
Q B 
Fia. 6. 


the upper half of the flow pattern is shown. As in Sec. 1, momentum considerations give 
the drag coefficient. If we normalize as before so that v, = 1, a = w/v, the drag coefficient 
based on the upstream velocity becomes 


“e = ] ( l — 2v ( ‘Os a + v’) y- / 23) 
C; = ; (i = 2v Cos a + v’) v. (23’) 
} 


Using (9.2), we get almost immediately 


dz = dz, + dz, 17b) 


where dz, is given by (17) and dz, by (17a). We therefore have b = b, + b, , where b, 


ta 
is given by (21) and b, is given by (21a). Moreover in general 
+ vf l o+f 1 + 6 Ll, 1+ %/8 


— =f log —— — 8 log — —- — log — —— 
—ve v v—¢ "1—p6F B °1—¢/8 


l 
z =v log = 
of 5 


where the origin of z is at the stagnation point on the lamina. Consider the channel 
wall CC’; along this wall ¢ is real and may be replaced by u, v S u S 1, and the complex 


coordinate along CC’ is 


; 1 + vu l v+tu iw ; 
z2=>2+i07 = 29 log + log - + — (A + A®). 
—_ ee l fe = FF i 
where 
Il+eu 
A =e” log 
l—e-u 
and A= is the conjugate of A. In particular at C’, u = 1 so that the distance from the 
plate to the channel opening is 
l+v, ] l+yv 
z. = 1 =v log 7. log ce tg (A + A*). 
== 7 / rr Fs 
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When u = 1, one finds that 


A + A* = 2 cos a log cot a/2 — rsina 





























and 
l= (» + 1) log x : + rsina — 2 cosa log cot a/2, 
‘ (23’) 
' = 2(: + 1) tan hv + rsina — 2 cosa log cot a/2. 
$= plane 
¢ le yn 
oO’ = l- yn plane 
it ¢ B,C" 
S O B,C C' 
Ss 
Fia. 8. 
T - plane 
T =-t2 T +0 T*1 T =s* 
—_ 
Cc’ c;B" ) .@) B,C C' 
; 
yrT yr yr" 
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A positive value of / corresponds to a plate position in the region of the free jet, a negative 
value to a plate position in the region of the closed channel. 

Graphical evaluation of the relation between a and v as given by (23’’), has been 
carried out for 1/(2r/v) = 0, + 1, + 2; and the correspondence of C, and C, on b/(x/v) 
for each of those values of 1 is shown in Figs. 10a to 10d. The dependence of the angle, a, 
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of the jet deflection on v is shown in Fig. 11; the dependence of the ratio, b/a, of plate 
half width to channel half width on v is shown in Fig. 12. Table 4 summarizes the nu- 
merical results. 

It is interesting to see what could be determined from tables of the incomplete beta 
function, in the cases of a wedge in free and bounded jets. By (20a), in a free jet, one 
could obtain the dependence of the wedge size b on the jet angle a, since | 8| = 1 and 

¢! = 1 at the points where the free streamline begins, knowing (18) only along a con- 
tour consisting of the real axis and the unit circle. By (17b), the same procedure is 
applicable to wedges in bounded jets, and one can locate the orifice as well. 

However, it seems impossible to obtain the pressure distribution along the tunnel 
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wall and dividing streamlines in these cases, without tabulating (18) over the entire 
complex plane. Hence we conclude that, with general n (not if n = 2), the fixed wall 
case is substantially easier than the free jet case, which is about as hard as the bounded 


jet case. 

9. Stability of the pressure coefficient. We can sharpen the conclusion of Sec. 7. In 
the fixed wall case, for any n, the theoretical pressure coefficient is approximately inde- 
pendent of v, if based on the velocity at the separation point, (downstream velocity). This 
we call the principle of stability of the pressure coefficient. 

TuHeorREM 5. In the fixed wall case, the pressure coefficient’* along the wedge, has a 
= |, for given wedge dimensions. 


zero derivative with respect to v, at v 
Proor. By (7) and (9.1), in the case of a wedge in an infinite stream, the complex 


potential 7 satisfies 


r= (14+ 2)’/1 — 2"). (24) 


18The pressure coefficient is 2(p — p-)/pv*, where p, is an ambient (cavity) pressure, and v an ambient 
velocity. Our proposal that one should use the free streamline velocity for v, instead of the free stream velocity. 
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Hence if we write v = 1 — e, we have 
r(v) = (1 + v")?/(1 — vo")? = 4/n?@(1 + ---) = O(1/e). (25) 


Now consider the functional relation between 7 and the complex potential W for 
the fixed wall case, for the same ¢. (They have the same hodograph; hence this is schlicht.) 
We normalize so that (in Fig. 1)W(S) = 0 and W(0) = 1—-e., so as to keep the total 
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change in W along the wedge constant. This makes the W-graph a strip of unknown 
width az, whence by (8) 
W = aLn (Ar + B)/(C7 + D). (26) 


We solve for a, A, B, C, D by matching W and 7 at ¢ = 0, 1, », e**””. This gives explicitly, 


writing 7(v) = 2, 


W = In (- = ) / tn (- a 


(27) 
if ( or) f ot /Qx* cee 
Ai + (2/22) + (2°/82") +) Le gy), 
\l + (1/2x) + (1/32°) + --: 
From (27) we conclude that changes in z = f dw/f are 0(e€) = O(1 — v)’ for corre- 


sponding values of ¢ and hence p, completing the proof. 
This suggests that in practice, cavity Cp should be computed on velocity at the 
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separation point, and not on upstream velocity. In estimating velocity at the separation 
point, we can either (i) use the Bernoulli equation and measure the cavity pressure, or 
(ii) assume uniform velocity in the section of maximum cavity diameter, and estimate 
the latter’s cross-section.—This discussion is only applicable physically to the case of a 


clear cavity. 
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THE CAUCHY RELATIONS IN A MOLECULAR THEORY OF ELASTICITY* 


BY 
IVAR STAKGOLD 


Cruft Laboratory, Harvard University 


1. Introduction and summary. If an elastic body is considered as a continuum, the 
existence of a “strain energy density” can be ascertained. 

Using a molecular theory of elasticity, we can deduce the existence of an energy 
function of the internal forces between particles. It will be shown that this energy 
function can be expressed as a function of an adequately defined strain. In the case of 
a monatomic crystalline body, this energy function corresponds to the strain energy 
density of a continuous body. For a multiatomic crystalline body, an analogy with the 
strain energy of a continuous body can be drawn. 

But in all these cases, it will finally be possible to write the energy density as: 


@ = $(¢;,). 


We can then assume the existence of a suitably convergent series expansion for ¢(e;;), 


d= bd + O'e;; + A” Ci iCnk + °°", 
(LD 


=g% +o tot: 
(see [1]t p. 230). The usual summation convention for an index appearing once covariantly 
_ ‘ — - ij. Ak . ° . 
and once contravariantly is used here. The coefficients A of the quadratic term in 
the expansion are those involved in the so-called “Cauchy relations.”’ 


We observe that A’’’ is a three-dimensional tensor of the fourth rank, and, as such, 
can have up to 81 independent components. By inspection of (1.1), we notice a number 
of restrictions on A‘’’“ which immediately reduce the number of independent com- 


ponents. 
Let us call (7, 7) the indices of the first group and (h, k) the indices of the second 
° > ij.Ak ° ° 
group. It is clear that we can choose the coefficients A‘’'”” to be symmetric in the two 
groups, 1.e., 


Atiem ai gee: (1.2) 


” , ‘ ia a ‘ vhk 
Since e,; is a symmetric tensor, it is also evident that we can choose A'’’” to be sym- 
metric with respect to indices of both the first and second group. Hence: 


gaveee = arrose - | deat = Da mess. (1.3) 


Taking into account these additional conditions on A, we conclude that only 21 inde- 
pendent components remain from the original 81. 
The Cauchy relations are not included in (1.2) nor in (1.3), but congist of six addi- 


tional relations 


ere = geniee ie preene ar (1 4) 
*Received Jan. 25, 1949. 


tNumbers in the brackets refer to the bibliography at the end of the paper. 
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whose validity we are to investigate. In other words, we are to find out whether A‘’*” 
is completely symmetric or only partially symmetric to the extent given by (1.2) and 
(1.3). 

This formulation of the problem was given by Brillouin ((1], p. 234) and previously 
suggested by Voigt ([5], p. 607) and Born ((2], p. 554). 

Within the framework of a continuum theory, no relations of the Cauchy type can 
be arrived at. It is only when we picture the elastic body as a set of particles with in- 
teracting forces that we can hope to arrive at the Cauchy relations. Naturally, we expect 
these relations to hold only under certain specific assumptions as to the type of body 
and the type of interacting forces. 

It is the purpose of this paper to discover the assumptions necessary (or sufficient) 
to insure the validity of the Cauchy relations. 

2. Historical summary of work on the Cauchy relations. As has been pointed out 
previously, all the investigations on the Cauchy relations must be based on a molecular 
theory of elasticity. 

Born was the first (at least to the author’s knowledge) to discover the fundamental 
distinction between the monatomic and the multiatomic cases. All previous authors 
had assumed that the concept of strain, as defined for a continuum, could be applied 
without any modifications to a molecular body. Born proves conclusively that this 
transfer is permissible only for a monatomic body, where each particle is a center of 
symmetry. He shows that the concept of a continuous strain could not be used directly 
for the multiatomic body. The other authors discussed in this chapter do not make this 
distinction, and their work is, therefore, only applicable to monatomic bodies. 

The Cauchy relations were claimed by Cauchy to be valid for the case of a molecular 
body in which the interactions between particles are central forces. Many writers have 
verified the Cauchy relations in the case of central forces and have shown that these 
relations are not valid if the interaction forces are not central (see [2], p. 555; [3], pp. 
47-49; [5], pp. 607, 608). All these verifications (except Poincaré’s) are based on a linear 
approximation to the strain tensor, rather than the correct expression for the strain 
tensor. The error, thereby introduced, will be shown to be of the order of magnitude of 
the remaining terms. 

The method of attack used by the majority of authors (including Voigt and Born) 
involves a series expansion for the energy function ¢ in terms of the Cartesian displace- 
ment components of the individual particles. After such an expansion is obtained, the 
displacement is expressed in terms of the local strain and is then substituted in the 
expansion. The resulting series for ¢ is a development in terms of the strain. The co- 
efficient of the quadratic term is examined to ascertain the fulfillment of the Cauchy 
relations. 

The difficulty of the method lies in expressing the displacement in terms of the strain. 
This involves the solution of the set of non-linear partial differential equations of the 


first order, 


1 fds Os; ' ] OS, OS, 

=5(=+—)+ Se Se (2.1) 
9 ee 4 . @ Te a 
«= \oy OY / oh CY CY 


Within the region of appreciable molecular interaction, it can be assumed that e,; is 


constant. 


1950 CAUCHY RELATIONS IN A MOLECULAR THEORY OF ELASTICITY 171 


In the one-dimensional case, (2.1) can readily be solved: 
s = (1 + 2e)'’x — z, (2.2) 


which reduces to s = ex for very small strain e. It is, therefore, possible that a closed 
solution to (2.1) might exist. The exact solution must have appeared too laborious to 
writers on the topic, for no recorded attempt at an exact solution can be found in the 
literature. : 

3orn and Voigt neglect the quadratic term in the derivatives in (2.1) and arrive at 


8; = ey". (2.3) 


Expression (2.3) is substituted in the series for the energy density ¢. Born and Voigt 
then conclude that, in the case of zero initial external forces, the Cauchy relations are 
satisfied by the coefficient of the quadratic term in the strains. 

Epstein, in 1946 (see [4]) presented a justified objection to the approximation (2.3) 
used by Born. He states that (2.1) must be solved up to terms quadratic in e;, to insure 
derivation of the complete quadratic term in the strain series expansion for ¢. He points 
out that this more delicate approximation need only be substituted in V, , the linear term 
in the displacement series for ¢. In the quadratic term in the displacement series, the Born 
formula (2.3) may be substituted. 

There can be little question as to the propriety of Epstein’s observation. Unfortu- 
nately, he uses an unnecessarily complicated method to arrive at the conclusion that 
the Cauchy relations are satisfied by the coefficient of the complete quadratic term in 
the strain series for ¢. Then Epstein states that the “‘coefficients of elasticity” do not 
satisfy the Cauchy relations. He considers the case where the initial external forces 
are zero. In this case, V; = 0, hence the energy is given by the term quadratic in the 
displacement (V,). He then states that this term is not equal to the entire quadratic 
term EF, of the strain series, and that, in addition, V, does not satisfy the Cauchy rela- 
tions 

In the case of zero initial external forces, the only case which Epstein considers, we 
will show that 2, and V, are equal and that the coefficients appearing in them are also 
equal. We will also show that, from a theoretical point of view, it is much more logical 
to define the coefficients of elasticity from 2, than from V, . Epstein does not notice 
that 2, and V, are equal. As a matter of fact, he explicitly states that his entire argument 
rests on the fact that they are not equal.* 

To arrive at our results, we will use a method of attack different from that of the 
majority of writers on the subject. This method has several advantages: it lends itself 
to simpler and more rigorous treatment; it does not involve approximations to the 
solution of partial differential equations; and it presents a clear tensorial outlook. The 
tensorial characteristics can be recognized immediately and all invariant relationships 
become apparent. These advantages are achieved by using, almost from the start, a 
series for ¢ in terms of the strain rather than a series in terms of the displacement. 

We will resort only to approximations of a physical nature, based on universally- 

*This has been pointed out by Zener (Phys. Rev. 71, 323 (1947)). A somewhat more detailed analysis 
is given by Per Olov Léwdin (A theoretical investigation into some properties of ionic crystals, Upsala, 1948). 
This latter pamphlet discusses the Cauchy relations in the case of finite strains and arrives at conclusions 
entirely in agreement with those in Sec. 4. The pamphlet reached the United States in the fall of 1948, 


after this thesis had been completed. 
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recognized experimental facts. These approximations are part of all other treatments 
of this subject, although the various authors have not always stated them clearly. 

3. Asummary of a few well-known results of the theory of elasticity for a continuum. 
Let the original coordinates of a point (i.e. before deformation) be represented by a’. 
The final or strained coordinates will be denoted by 2'. The coordinate systems need 
not be the same. 

Let ds, denote the square of the unstrained distance between two neighboring points; 
let ds” denote the square of the distance between the same two points after deformation. 
We then introduce the Lagrangean strain tensor e;; for the deformation of a continuous 


body; from its definition 


ds*° — ds, = 2 ij da‘ da’. (3.1) 


The existence of a strain energy representing the potential energy of the deformation 
has been shown by many writers. (e.g., see [1], p. 229). 
Let E be the total strain energy for a continuous body. We write 


E= | ¢dav, 
“1 
where ¢ is the strain energy density and dV is an element of volume in the unstrained 
coordinates. 
Since the strain e;; completely determines the energy of the body, it is clear that 


@ = ¢(e;;). We assume that 
2 es to, tg tees =H te; +A eeu tee. 3.2) 


Now ¢ is a scalar, as are go , ¢; , ¢2, etc. Hence 2"’ is a doubly contravariant symmetric 
tensor, and A‘’’”’ is a four times contravariant tensor symmetric in (7, j) and in (h, k). 
The two groups (77), (hk) are also interchangeable. No additional symmetry in this 
tensor can be discovered by treating the body as a continuum. Hence the Cauchy rela- 
tions are not satisfied. 

4. The monatomic lattice. Consider a simple monatomic lattice. All the particles are 
similar and have the same force law between them. We will assume, for the present, that 
knowledge of the position of every particle in the body completely determines the energy 
of the body. We are excluding the possibility of dipole moments attached to the particles; 
in this case, the energy would also depend on the angle between the dipole vector and 
some reference frame. It may seem that we are also excluding the monatomic lattice 
in which the particles are considered as rigid bodies with ability to rotate. It was shown 
by Born that this case reduces to the case of a multiatomic lattice (see [2], p. 556), 
which we will treat in Sec. 5: 

We conclude then that the total internal energy of the body depends only on the 
distances between atoms. This is not equivalent to the central force assumption, where 
the total energy is the sum of mutual energies between pairs of particles, this mutual 
mutual energy being a function only of the distance between the two particles on which 
it acts. 

Let u, v, 7, *:- be particles of the body, p,, the distance between the particles yu 
and », £ the energy function for the internal forces. It will be more convenient for our 


ge 
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purposes to use R,, = p,» aS the independent variable; this is permissible since p,, is 
positive. 

E = E(R,,, Ry, Ry» , +++), in general. (4.1) 


The case of central forces occurs when the function E can be decomposed into 
] . , 
»s > E,(R,.), for central forces. (4.2) 


(In (4.2) R,, , R,, are treated as independent variables.) 
Expanding (4.1) in a power series, we get 


E => EAR,» ; Cn ’ “9 + > (22) are 


ue 


(4.3) 





OE -) 
+5. p> (a2 OR: AR AR, + - 


win.€ 
where the summation is over all pairs of particles 4, ». We may take E, = 0 as our 


datum. 





Let us divide the body into regions V, V’, --- . We define the energy of V as 
7 OE ) (22 ) 
E, a, 2 (2 AR,,. ne. 2 2 : OR,, OR, AR AR, + ° (4.4) 
in V 
We notice that E # y, Ey, , but that we have 
Eo: > Ey. (4.5) 
all V 


It is shown in the appendix that the equality sign in (4.5) becomes valid as long as each 
volume V has dimensions which are large compared with the distance of appreciable 
molecular interaction. 

Let us now consider a continuum superposed on our monatomic lattice before 
deformation. Each particle yo will have a coinciding representation y,, on the continuum. 
The continuum upon deformation will have a strain e;; , which is a continuous function 
of the coordinates. Hence, within a small volume V (which can still contain a very large 
number of representations u,), the strain is sensibly homogeneous. 

For a monatomic lattice, equilibrium is achieved in the unstrained state by having 
each particle a center of symmetry (the laws of force between all particles are the same). 
Are we permitted to assume that, after deformation, the particles still coincide with 
their representations on the deformed continuum? The only way the assumption can 
fail is by causing violation of the equilibrium of individual particles after deformation. 
But under the proposed assumption, each particle » will still be a center of symmetry 
after deformation (at least with respect to the nearby particles, which are the only ones 
influencing the force on «). Hence, we conclude that the assumption is permissible for 
a monatomic lattice, since equilibrium of all particles is maintained after deformation. 

Consider a volume V whose dimensions are small macroscopically speaking, but 
very large compared with the distance of appreciable molecular interaction. (From our 
knowledge of the forces of interaction, such volumes exist.) The assumption of the 
existence of such a volume V about each particle of the body is necessary to treat 
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inhomogeneous strain (this assumption is made by all authors, although not always 
stated). 
Within JV, e,; is nearly constant. From the definition of €;; , we have 


AR >, = ey 0,0, : if MM, v are close together (4.6) 
(e;; 1s the Lagrangean strain tensor; aj, = a; — a, , where a; are the unstrained co- 
ordinates of the point u). Since V is a small volume on a macroscopic scale, (4.6) holds 


for any two particles of V. 
Substituting (4.6) into (4.4) we obtain 


/ ‘ OF \ l OE ‘ , 
Ey = } 2e;;a),a), + = 16; CnpQiyA iy Agent vee 1.7) 
2 \OR,,,/ _— ~ 2 \OR,, OR y/o purge Ang T 
n } in ¥ 
We are now in a position to define the strain energy density ¢ = limyio (Zy/V). The 


limit is taken macroscopically, so that the smallest volume still contains a very large 


number of particles: 


] [ ' OE / l GE | “a 
eo" Vv L 2» Poy | 26st + 5 2 aR ah.) nee | os 


i | 


Now (4.8) is a series of the type (1.1) and we immediately pick out 


2 ( @#E , 
A = — : ; -] A, A, yAntAne |, (4.9) 
J p> : \aR.. OR, ¢/ * SiN iad ia 
where (0°£/0R,,0R,;), is a scalar in the unstrained coordinates. In general, A will 


not satisfy the Cauchy relations, since indices of the first group are not interchangeable 
with indices of the second group. For instance, i and h cannot be interchanged in (4.9). 
If we assume central forces, then substituting in (4.7), we get: 


E= : 2. | a 


1.10) 
’ ‘dz... IE ; 
Br = ¥ (the) vata +S (Ge) cena b> 
for central forces) 
and 
hi ] ‘a°E i j h } 
ie = , >» (see) aaa ; (4.11) 


ie 3 
It is apparent that the Cauchy relations are satisfied, since complete symmetry in i, j, 
h, k exists. 

It should be noted that no approximation to the strain tensor is made throughout 
the proof. The strain tensor is introduced through its very definition. The only approxi- 
mation used is one which is used by all other writers, i.e., about each internal point of 
the body there exists a region whose dimensions are large compared with ranges of 
atomic interaction, but in which the strain can be treated as homogeneous. 








~I 
or 


1950) CAUCHY RELATIONS IN A MOLECULAR THEORY OF ELASTICITY 1 


It is also well to remark that the coordinate system used in the proof is not specified 
and need not be a Cartesian system. Hence it can be concluded that (4.9) and (4.11) 
are, indeed, general tensor formulas. 

Having arrived at the desired results, we are now in a position to compare our work 
with that of other authors on the subject. There have been so many authors who have 
written on the subject that it will be impossible to discuss the investigations of all of 
them. The classical proof of Cauchy’s relation is well given in Born ((2], pp. 545-555). 
We shall limit ourselves to a discussion of his work. In addition, we will look into the 
work of Epstein as presented in his article, which was one of the impetuses for the 
writing of this thesis. Since these authors restrict themselves to central forces, we will 


write 
: 1 /dE,, ’ @E,, : _ 
y= Dig (i) aR X (ae ) ar.) ite ais 
in V in V 


Equation (4.12) is equivalent to (4.4) for central forces, and yields (4.10), from which 
we deduced the Cauchy relations. 

To conform with the work of other authors, we introduce Cartesian coordinates y’, 
and the displacement s; . Distinction between covariant and contravariant components 
is now immaterial. 

Using Cartesian coordinates y}, , we have 


AR,, = > (yj, + 81)? — Dd (yi,)” 


2Y Sine + rH (Bias 


(4.13) 
- ZY Sino + 5°'8; Sine ° 
We shall use the relations 
08.) _ (ibe) oy 
= 0 - dR,» Use ; aad 
aE, ) “(ae ) F (zz ) 
Te) ae 25'S) + ayy! | Soe) | 4.15 
F OY» 0 = dR,, 0 * YurYus dR, 0 ( ) 
Substituting (4.13), (4.14), and (4.15) into (4.12), we get 
. 1 IE, l OE,» 
Ey = 9 Ze (22) Sine 4 +> (- * oy ) sai pte s 
in V OYur 7 in V Yur Yur . (4.16) 
=VitVt-, 


where summation is over all repeated indices, regardless of position. Equation (4.16) 
is exactly the formula used as a starting point by Born and Epstein. 

Now, Born and Epstein both start from (4.16) and try to reach (4.10). Born is un- 
able to do this correctly since he does not use the correct expression for s,,, in terms 
of €;; . 

Born substitutes in (4.16) the linear approximation s;,, = €4;Yj» : 


, 1 0E,.\ ; 1 vE,, . 
E, = 5 > (see) wise + 4 > (<7 Be) eaenytott wee, (4.17) 
eis in V 


, OY i» OY ur 
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Using (4.14) and (4.15), he arrives at 


: 1 dE,, bos l dE 4»\ o sii 
Ey = 5h (iB: ) eaters +42 (se) 29 YarYart xeon 
7 


in in V 


(4.18) 


l VE,, i h 
+ 4 p> (se ) duySvbusevtsewen sg 
in V 


By comparing (4.18) with (4.10), we see that Born’s first order term is correct, but that 
his expression for the second order term is incorrect by a factor 


1 (2 ) bk 

Uh! = 25°i{ =") ¢..e.,y'y*, (4.19 

. 1 2d dR,,/o Winter ) 
in J 


(see [2], p. 547, eqs. (31) and (32)). 
Notice that (4.19) does not satisfy the Cauchy relations from a tensorial point of 
view. How then is Born able to deduce the Cauchy relations? He notices that, 7f the 


initial external forces are zero, then Ei’ = 0. Born actually never states explicitly that 


E;/ = 0, when the initial external forces are zero, but he uses this result repeatedly to 
simplify his equations. Kellerman, a student of Born, does make the explicit statement 
that £3’ = 0 and attributes it to Born (see [8], p. 534, note). In this case expression 


(4.18) reduces to the correct expression (4.10) up to terms of the second order (terms 
of higher order will always be incorrect). This is really a very fortunate coincidence, 
* since it led Born to the correct conclusion. 

In the case of initial external forces different from zero, Born would have arrived at 
the incorrect conclusion that A‘’’” did not satisfy the Cauchy relations. Expressions 
(4.10) and (4.11) show that the Cauchy relations are satisfied even though the initial 
external forces be different from zero. 

Let us now examine the work of Epstein. He considers only homogeneous strain for 
a monatomic crystalline substance (although he does not state this explicitly). He starts 
from (4.16), which he first rewrites in an unnecessarily complicated manner. He in- 
tegrates approximately the partial differential equations connecting the strain and the 
displacements (1.9). His final formula for s;,, is: 


Sine = SisYou — FSi SinYar » (4.20) 
(see [4], eq. (14)) where f;; = e;; + ,; and w,; is the rotation 
= ] (281m deur) 
wt 2X ay, yi! 
(Born uses 
Sige = Cis» +) (4.21) 


Within the volume J, it is clear that e,;; , w;; are constant. He then substitutes (4.20) 


into V, in (4.16). In V, (the term quadratic in the displacements) he substitutes s;,, = 


SisYur , a8 Born does. 
Upon substitution of (4.20) into V, we will get a term linear in e,; , call it V{ , and a 
term quadratic in e;; , call it V{’. Epstein considers V{’ + V, (which is the whole term 
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quadratic in e;; , and which we have previously called g, (energy density) or E, (energy)). 
He sets the coefficient of the expansion of V with respect to the powers and products of 
w;; equal to zero. He then arrives at a number of rather complicated relations between 
the coefficients, which he uses to simplify the terms involving only e;; . He then con- 
cludes, after this laborious task, that the coefficients of [V{’ + V.] (=E,) satisfy rela- 
tions of the Cauchy type. This we have proved in a simpler and more straightforward 


way 

He restricts himself to the case of zero initial external forces, and shows that V, = 0, 
hence that 2 = V, + --- . If we neglect higher order terms, we have E = V,.(V,; = 0 
since 


ee 
vy \OY,,»/ 9 
is the ith component of the initial force on the particle u = 0.) 
He then has two quadratic forms in e,; at his disposal: 
(a) Vi’ + V., which satisfies the Cauchy relations (= F,), and 
b) V., which is the energy. 
He now claims that the coefficients of elasticity appear in V, and not in Vj’ + V,, 
hence that the coefficients of elasticity do not satisfy the Cauchy relations. Actually 
the coefficients of elasticity (for the case of zero initial external forces) appear both in 
Vi’ + V, and in V,, since it is easy to show that V{’ = 0. 
We have that 
1] ij7,hk 
@ = Q'e;; + A’ C5;€ne tee 
Since the original state is stable, it is clear that Q'’ = 0 for the case of zero external 
forces. This can be shown in another way, using a formula which Epstein takes as a 


basis for his argument: 


i= a Of + 2AM, + oes, (4.22) 


ov ij 
If the original external forces are zero, then the initial stress is certainly zero; hence 
Q'’ = 0. Consequently 


¢=¢d+°::: or E=E,+°:::. 


Limiting ourselves to terms of the second order in e;; , we have that 


E = E, = Vi’ + V,. (4.23) 
But we also have that EH = V, . The only possible conclusion then is that V7’ = 0. 
Formula (4.22) shows clearly that if any terminological distinction is to be made be- 
tween the coefficients in FE, and V, , the coefficients of FE, = V{’ + V., and not the 
coefficients of V, , must be called the coefficients of elasticity. In any event, the co- 
efficients are equal and E, = V{’ + V,, represents the energy just as well as V, . 


The fact that V{’ = 0 has been conclusively proven by Poincaré and Born (e.g. see 
[3], p. 45). It might be illuminating actually to substitute (4.20) and (4.21) in their 
respective positions into (4.16) and see that Epstein’s formula for the energy is indeed 
correct (of course only up to second order terms, which is all he claims). Actually, 
Epstein goes about this in a more complicated way but with no apparent advantage. 
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We have upon substitution: 


, l OE »\ - 21 ] JE, 
Ey =35 (2 Ja Yur — | >» (: “) Sirf rYor 


= OY)» ar OY; 





in V in wii 
(4.24) 
0°E, es = 
+ Oa : L -) SisFuniYurU re + nt es 
2 OY,» OY;,»/ 
Using (4.14) and (4.15), (4.24) becomes 
a (A) 1 (sz :) 2 
Ey = 9 a a dR,, 2YirYne. fis - 4 2 IR. YieYuod iat ir 
in V in 
| (dE : = 
= Qsti| ““*ur f 4™ 2! 
4 2 ” (SB Yurir SriSmi (4. ») 
in V 
d ft) 4 bot ee 2 
Tz iz iF | AY uYurYurYurdiiSmi 1 
or 
= l dE,, F _ 
Ev=5 2 (ee #) ytd 
in V 
(4.26) 


1 aE, 
$ 4 p> (E “) AY) YreYurY ert ti. iT Peas 


Setting all the coefficients of terms involving w,; in this equation equal to 0, we get 
(4.10) up to second order terms: 


ss. (iz. ‘) — 
E v™ @g p> dR,, 2YurYurl 


~ 


in v 


(4.27) 
| l (7B) ae ? h _k 
! 4 > dR: } FY rYurYurY wel i j@hk + ’ 
in V 


Ey =E£,+F8.+::: 


The quantity Z, Epstein calls V{’ + V, because, in his derivation, Z, stems partly 
from V, and partly from V, . This has no theoretical significance. We can derive (see 
derivation (4.12) to (4.16)) the series in the displacements from the series in terms of 
the strain and find that EZ, yields V, and part of V, . These statements are in no way 


contradictory. In any case, Fy = FE, = Vi’ + V. = Vz (see (4.23)), where the initial 
external forces are zero. We have the expression for V{’ in (4.25): 
: ] dE 
Vr’ = — = ( ) Dit Obs (4.28 
1 4 > dR,, é YurY uv Ke ik ) 


in V 
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The crux of Epstein’s argument is that V{’ # 0 (see [4], p. 919); he states that the only 
difference with Voigt is a term 


1E 2\2 3 \2 
> (see) ut, + (Ys»)'. (4.29) 
inv liad 


Epstein also says that this expression does not, in general, vanish. Yet it was proven 
by Poincaré ({3], p. 45) and by us that V{’ = 0. Hence (see (4.28)), we have: 


a) i i = 
»» (33 Yuan = 0 (4.30) 
in V 


(zero initial external forces). This also proves (4.19) used by Born. Hence certainly 
(4.29) = 0, and the entire argument of Epstein breaks down. 

We have Ey = E,(=V?’ + V,), and the coefficients of Z, are the coefficients of 
elasticity and satisfy the Cauchy relations. 

5. The multiatomic case. Let us consider a multiatomic, crystalline, elastic body. 
It can be thought of as a superposition of many simple monatomic lattices. The body 
can also be represented as the repetition of some fundamental cell (containing all the 
different types of particles!). 

It has been shown by Born ((2], [6]) that, if the laws of force between the particles 
of different simple lattices are not the same, it is impossible to consider the body as 
strained in the same manner as a continuum. He shows that this latter assumption would 
lead to a violation of the equilibrium of the particles (see also Love [7], note B, pp. 
620-627). 

We can make a very reasonable assumption as to the state of deformation of the 
body without violating the equilibrium conditions and without abandoning the concept 
of strain as defined for a continuum. Each individual monatomic lattice is strained as 
a continuum, in the manner described in Sec. 4. In addition, each simple lattice is dis- 
placed with respect to the others. 

Let e;; be the strain of all lattices. Let us picture a frame of reference consisting of a 
continuum C, superimposed on the molecular body, and let this continuum be subject 
to the strain e;; experienced by each of the monatomic lattices. Every particle u of the 
body will have an initially coinciding representation yu, on C. The representation of 
each particle on C will be displaced under the strain e;; by an amount s,, . The final 
position of the representation yu, of a particle will not, in general, coincide with the actual 
final position of the particle u, since each lattice is undergoing an additional displace- 
ment with respect to the continuum. 

Let the displacement between the deformed positions of a particle and its repre- 
sentation on the continuum be wu, . In a region which is small on a macroscopic scale, 
ii, will be the same for all particles « belonging to the same monatomic lattice. If our 
body consists of n different lattices, the deformation is completely specified by the 3n + 6 
local constants 


ua es U2): ia SS U(m)i “Wixthging'y Un) i , Ci; ° 


We recall that the u,,.; were introduced to insure preservation of equilibrium. We 
expect the imposition of the equilibrium conditions to yield u,,,); as a function of the 


local strain e,; . 
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In any event, let us start our development as in Sec. 4 and notice the differences as 
we go along. For a monatomic lattice we expected an energy density of type (1.1), 
which was easily predictable. For a multiatomic body, we expect the energy function 
to depend at first on e;; , Ua), *** » Us - We write tentatively 


hk 


0765; FA ene + 
2 WimyUimys + : ms ny cm? )Uny tmy3 


g¢=>o 7 


+ 
(5.1) 
+ DT istewmess + 
We have that 
Od 
n= [58] <0 
ou m . 0 
to satisfy the original equilibrium. For final equilibrium, we must have 0¢/d0u;,); = 0’ 
which leads to: 
SS. Bw tian + Tate @ 4. (5.2) 


We have 3n linear equations in the 3n unknowns vw,,,); . These equations are not inde- 
pendent but still admit of a solution (see Born [2], p. 552). The solution will yield 
as a linear function of all the strain components. We expect that 


Um) ij 
Wim), = Winlhk « (5.3) 

Substituting in (5.2) leads to 
o> Buwtesta + Tits = @. (5.4) 


Using (5.4) we simplify (5.1) to give 


1 Ak 
@ = Q'e;; +A eCijn + eee 


se H m 1... jOnxlip Te? 


The terms appearing in (5.5) will be evaluated a little later on. We will notice then that 
A‘’’™ satisfy the Cauchy relations under the assumption of central forces. We will also 
observe that the first line in (5.5) is identical with the energy density based on an assump- 


tion of continuous strain 
The complete coefficient of the quadratic term in e,; will satisfy the Cauchy relations 


if, and only if, 


i eS AE A ro (5.6) 


satisfies the Cauchy relations. 


, ‘ . . . tj,hk 
We recall that the Cauchy relations held in the monatomic case because A 


could be written as 


~] 


et 


A = ‘> A,,@;,,0),0,,Q,, (A,, = scalar). (5. 


uy? 





1950) CAUCHY RELATIONS IN A MOLECULAR THEORY OF ELASTICITY 181 


It is clear that the simplest tensor of the fourth rank satisfying the Cauchy relations is 
xi = v'v'y'y* (5.8) 


where v* is an arbitrary contravariant vector. 

Let us try to make the most optimistic assumptions possible in order to obtain the 
Cauchy relations. Let us assume the form of 77;,,; most likely to yield a B“’’” satisfying 
the Cauchy relations. Clearly the most favorable assumption is 


hk h k 
Timi = Us (my b(n) b(m) > (5.9) 


where ¢?,,) is a vector depending on m, i.e. different for different m. (Notice that x°*,, is 
a quantity related to the mth lattice, and must be assumed different for different m.) 
Let us assume for simplicity that 


‘7 i i - ~ 
Hn) (m") = Qim)yQ@(m*) ’ (5.10) 
then 
hk, lt ? h k I p - 
B = > (my (an?) bi (my bi (my EG my Em) Ean") (m') (5.11) 


»(m") 


In addition let us assume that the summation over (7, 7) can be performed, to yield 


BS? alone 2 Colton (5.12) 
(m),(m") 

By inspection of (5.12) we see that the Cauchy relations are not satisfied as soon as we 
have two or more types of particles. For, let us write 

Bp" «= (ee) 2. GishiaSeiheas « (5.13) 
Then comparing (5.12) with (5.13) we see that the coefficients are not equal since 
tim) * time) at least for some values of (m), (m’). Hence the Cauchy relations must 
fail for a multiatomic body even though subject to central forces. 

It may be of interest to develop in detail the expression for ¢ in the case of a multi- 
atomic body, thereby verifying, step-by-step, the argument given previously in general 
terms. We then can compare our series for ¢ (or E) with the one used by Born. We will 
arrive at the conclusion that the series used by Born is incorrect in the part involving 
the strain terms exclusively, but correct in the terms in u;,,,, and the mixed terms in 
Uscm)Cnk » 

We will limit ourselves to the case of central forces, and we start with formula (4.4), 


rewritten for central forces: 


2 28 dE,, i (<Ee) Dis on ‘ 
Ey =} > (B+) an + 4D (Get) (aR. + oo (5.14) 
in V in V 


Unfortunately we cannot write, as we did before, 
AR,,, = 2€;;QprAye - (5.15) 


Formula (5.15) holds only if u, v belong to the same monatomic lattice. If u and v belong 
to different lattices, then, in addition to the strain e;; , we must take into account the 


difference in the additional displacements u, and %, . 
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Let s, be the actual displacement of the particle u. Let s,, be the displacement of 
the representation of u on the continuum C. Then 


Ss = §. + Uu,, (5.16) 
and 
8.» = 8,.¥. T Us » where ti,, = ui, — u,. (5.17) 
Let R,.,., = Pusr. » Where p,,,, is the initial distance between yu and p; 
Ru.r. = Pucre » Where p,.,, is the strained distance between the representations of 
uw and y on the continuum C; and 
R,, =p» where p,, is the deformed distance between wu and p. 
p 
The following diagram may clarify the situation. 
4 
Oo 
~ 
~~ 4) 
~~ © 
Sy vy, \ a 
Ko \ ~ v, 
in 
| \= pene [ 
“ Yn Fhe%e 3) | 
/ He \ 
| & \ | 
Lae 
/ V] 
Pr Puv 
We have 
AR,. = Ry. — Rare, » = Rue — Ryso. + Rucve — Reco ; 
whence 
Reeve ~ Reuscre = (AR) - = BWeijAirAus « (5.18) 
Clearly 
Kk, = &,.. = 2,2) + 2. (1,3); (5.19) 
where we are now using Cartesian coordinates. 
Consequently, 
AR,, = (AR,,)- + 8 Uy Ua; + Ys, + Sir.) Uyrs (5.20) 


(since Cartesian coordinates are used, s' = s; , etc.), or 


AR, = 2€:;YsYoe HO Uyy Ua; + ivi + yo ore (5.21) 


| 
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Substituting in (5.16), we obtain 


uy 
n b in V 


7 ] dE,, 9 (fe. “) 4 ' i a a 
Ey = 9 > (2) 20 ;Y, iiles a >> dk, ‘ $6 CY ur YurYurYur + : 


| dE \ IE UE _ 
3.3 (te) rina +3 % (Me) nas #3 & (BE) ata 
n V in V in V 


(5.22) 


pv 70 


[ss -) se 
oto QR c F 
) i > dR. / Luk YuvYurY uv a 


7. (tHe) ee 
T? » AR yy) 9 reretons TF 


In (5.22) we have included all terms up to and including second order terms in e,;; and u,,; . 

The first line of (5.22) has the familiar appearance of the strain energy encountered 
in Sec. 4 (see (4.10)). It represents the energy resulting from the deformation of all the 
lattices, excluding the energy due to the displacements of the lattices with respect to 
each other. Using an obvious extension of the terminology already introduced, we will 
call the energy terms in the first line of (5.22) by EF, . 

Examining the other terms appearing in (5.22), we first notice, by investigating the 
condition fot equilibrium, that 


a (aes Yutp»s = 0. (5.23) 


“PY uy 0 
in V 


We are next interested in the cross-terms involving e;;u,,, . The only direct term of 
this type appears on the third line of (5.22). A veiled term of this type, however, appears 
in the fourth line, which includes all cross-terms in s’,.,.u,,; . We know that s).,. = 
f(e;;), and since we are interested only in the term linear in e,;; , we can safely substitute 


] ik ? 
ie Cri Yur = 6 Cei Yur 


We are now in a position to write an expression equivalent to (5.1). We have not 
yet written u,,; in terms of u,,); , but this is only a simple algebraic step. We have 


j I ‘dE, — l aE i ee 
9 » (th) 9 i YurYur + | z (se £) des seused ett + 


| » J ‘3 v 
41] ss ib.) is (TE) | = 
@ = V ) + 2 > E (se 5 + SY uYur 1R? u prillyy; + | (5.24) 
in V 





| by ie) (E-:) | 
+ p> E yao s + 2 Qy. YY iR?, s Ci Upvk oe | 


in ¥ 
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By using relations (4.14), (4.15), we rewrite (5.24) as: 





be a Peary ” | 
= . = i“ Upvityr; 
| 4 a \ayi, dyi,/> 
} in V 
, — 
you eds f (5.25) 
1 oE,, 
;| +2> (: a) Ye 
mH wy OY uy OY u» 0 
in V 


We can arrive at (5.24) or (5.25) in a somewhat easier way by using the results for 
a monatomic body, and starting from a series in terms of the displacements: 


9B a 
. - k, + : - - | a l > ( ¢ Eu) St os, (5.26) 


OY.» i OY,» OY.» 


in } in 





Me? uy 
V 


With E, = 0, ands, = $,. + u,, 


‘] dE ,, l oE,, 
org | rer eeweg | ee 
vy \OY,»/ 0 t uy? OY,» OY.»/ 0 
J 


bo 


in in 


uv \OYL, OYi»/ 0 


.. 3 I (OE, ] or... ~ OF 
= V ao 5 2. (: = Unvs ~ 4 > oe) Use ithess a ae (5.27) 
OE,, \ 


+-> a( =e Sereitlues 
Bex OY pn» OY ur 
a J 


Now (5.27) is equivalent to Born’s series (see [2], p. 547, eq. (31)). 

We notice that the first line of (5.27) is exactly of the type found in the monatomic 
case and is equivalent to what we have called ¢, = E,./V. In the third line of (5.27), 
we are allowed to substitute 


We notice that 


(22. 
Ua»; = 0 


wee OY )y/ 
in }¥ 
and that (5.27) therefore becomes (5.25). 

30rn substitutes s,.,.5 = €:;Y,» throughout (5.27); actually this substitution is only 
permissible in the third line of (5.27). Hence Born arrives at an incorrect expression 
for Ey and ¢, , but, since his error is localized in ¢, , he is still able to deduce correctly 
the equations equivalent to (5.2) through (5.6), which are the salient equations of the 
multiatomic case. 
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APPENDIX 


We have written the total strain energy of a molecular body as 


ina =F.) ( 7 oe | 
B= 2 iz Aor + oe NSB, alta) amare + ++» 
ins ' (A.1) 
=F,+EF.+-:-, 
where ; ae means a summation over all pairs of particles y, v. 
We then defined the strain energy of a volume V as: 
; a aE 
Ey = (2 ) arn + > (22) AR,,ARy + °°: , 
y \Ofbyy/o uivin, & \OLbyy Olbye/o 
“ V in V ' (A.2) 
Ey, + Ey, + 
If the body is divided into two parts V and V’, we notice that 
Ey a Ey: x E. (A.3) 


The reason for this inequality is physically clear. We have omitted the surface energy, 
i.e. the energy between the volumes V and V’. This energy is denoted by Eyy- . It is 
our contention that yy. is negligible compared with Ey and Ey. if the volumes V, V’ are 
large compared with the sphere of molecular activity. Then the inequality (A.3) would 
become an equality: 


Ey + Ey = E. (A.4) 


We will first show that Fy, + Ey, = FE, , the arguments for higher order terms running 
similarly. We have exactly: 
Evy + Ev + Evy = Ei. (A.5) 
We notice that Zy;.y-, consists of a summation over pairs of particles, one of which 
is situated in V, the other in V’. The only pairs Of particles, which will contribute appre- 
ciably to this summation, are those within an extremely narrow zone about the boundary 
of V and V’. The volume of this zone is of the order of the product of boundary surface 
multiplied by the range of appreciable molecular activity. This volume is therefore 
very small compared with either V or V’. Hence, it is permissible to write Ey, + Ey, = 
E, . The reasoning for 2, , etc. --+ is similar, and: 


Ey + Ey = E. 


If the body is split in more than two parts, we will still have 


se 
The above argument is given in great detail by Poincaré ([3], pp. 40-53), and we have 
done little more than reproduce his exposition. (Note that long-range forces are not 


included in this treatment). 
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CONDUCTION OF HEAT IN COMPOSITE SLABS* 


BY 
J. C. JAEGER 


University of Tasmania 


1. Introduction. Among the most important practical problems in Conduction of Heat 
are those of flow of heat through composite slabs of several layers with different thermal 
properties such as occur in the walls of furnaces, rooms, and refrigerated chambers. The 
well known result for steady flow of heat in such a region is very simple: if 1, , k, , and 
R, = I,/k, are the thickness, thermal conductivity, and thermal resistance, respectively, 
of the r-th layer, the rate of flow of heat, per unit time, per unit area, per unit tem- 
perature difference, through a wall of n such layers with perfect thermal contact between 
them is 


{R, + R, +--+ + R,}7', (1) 


and the same result holds if some of the layers consist of air spaces or surface films of 
negligible heat capacity provided each film is represented by an appropriate thermal 
resistance. 
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Fig. 1. 


In contrast to the simplicity of (1) the complete solutions of transient problems on 
composite walls,’ even in the simplest case of two layers only, are so very complicated 
as to be practically useless. Information is needed about many transient problems such 
as the warming up of a furnace, the initial cooling of a refrigerated chamber, and the 


*Received April 19, 1949. 
For references and some solutions see Carslaw and Jaeger, Conduction of heat in solids, Oxford, 1947, 


Art. 119. 
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thermal behaviour of buildings. Broadly, what is required may be either (7) a prediction 
of the transient behaviour of a given wall under given conditions, or (77) a comparison 
of the transient behaviours of two different walls. 

The first object of this paper is to discuss, in detail and from both the above points 
of view, one of the commonest transient problems, namely, that of finding the total 
quantity of heat which passes through unit area of a wall from zero time up to time f, 
the wall being initially at zero temperature and its surfaces being maintained at tem- 
peratures V (constant), and zero, for times t > 0. As remarked above, the complete 
calculation of this quantity would be extremely laborious, nevertheless it will be found 
possible by a simple arithmetical calculation to determine accurately a quantity char- 
acteristic of the wall, which we shall call the time-lag in establishing the steady state, 
which gives sufficient information for many purposes. 

To illustrate what will be done, we consider first the case of a wall consisting of a 


single layer of thickness /, conductivity k, density p, specific heat ce, and diffusivity 


a k/pc. Suppose that the wall is initially at zero temperature, and that for t > 0 
its faces x 0 and a / are maintained at constant temperatures, V and zero, re- 
spectively. Then Q) and Q, , the quantities of heat which cross unit area of the faces 
x 0 and x l, respectively, up to time ¢, are given by” 
a | ‘ 2kVl oo 1 ya 5 
@) i+: " ; exp (—anr t/l), 2) 
l 5a ar d i 
RY ( r) 2kV1l S (-1)" = 
) , == s—- exp (—an'ar' t/l). (3) 
; l 6a ar > n ditties 


For large values of the time the exponential terms in (2) and (3) tend to zero, and 
(), and Q, tend to straight lines of slope kV/l with intercepts on the time-axis of (—I°/3a) 
and (l°/6a) respectively. The way in which these linear asymptotes are approached is 
shown in Fig. 1 in which Q,/Vipe and Q,/Vipe are plotted against the dimensionless 
quantity at/l. It appears that the curves approach their asymptotes rapidly, and that 
at times of the order of the greater intercept, namely [°/3a, the departure from the 
asvymptotes is quite small. Thus for times greater than this, the values of Q, and Q, 
are the values for steady flow for times t — 7. and t — rt, , where 7, °/3a and 
T l’/6a may be described as the time-lags in establishing the steady state 

For the composite wall of » layers, the complete solution for the quantity of heat 


crossing either surface up to time ¢ will have the form 
Q) ViR, + --- R,}'(t — r) + negative exponentials, (4) 


and for large values of the time will be a straight line of slope given by (1) with intercept 
7 on the time axis. In Secs. 2 to 4 formulae will be derived for calculating the time-lag 
7 and it will be shown that the graph of Q tends to its asymptote in much the same way 
as those of Fig. 1 so that a qualitative idea of the behaviour of Q for times less than + 
can be obtained 

The same method can be applied to any other problem in which the quantity under 
discussion ultimately increases linearly: solutions of a number of such problems are 
given in Sec. 5, they include cases of constant flux at a surface, linearly increasing tem- 


*Barrer, Trans. Faraday Soc., 35, 628 (1939). Carslaw and Jaeger, loc. cit., Art. 117. 
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perature at a surface, ete. In all these cases a time-lag appears which can be expressed 
in terms of the parameters defined in Secs. 3, 4. The behaviour of an unsymmetrical 
wall for opposite directions of flow is discussed in Sec. 6. 

In Sec. 7 the formulae derived in the earlier sections are collected for reference and 
a typical calculation set out in detail. It may be remarked here that the formulae are 
expressed in terms of the thermal resistances and capacities of the layers, so that they 
may be applied immediately to walls containing air gaps or films of finite thermal re- 
sistance and zero heat capacity. 

The time-lags so obtained, being accurately calculable quantities which appear 
naturally in the solution of transient problems on composite walls, are also available 
for the comparison of the transient behaviour of different walls: it is suggested in Sec. 7 
that the time lag in the quantity of heat passing through a wall when the temperature 
of its outer surface is suddenly changed gives a satisfactory criterion for such comparison. 

2. The composite slab of n layers with perfect thermal contact between them. Surface 
temperatures constant. Suppose the r-th layer has thickness 1, , conductivity k, , specific 
heat c, , density p, , and diffusivity a, = k,/p,c, . We suppose the layers to be the regions 
0<a<l1,,l,<a2<1,+1,---, ete., and in this section and section 4 that the outside 
surfaces, x = Oanda = 1, + --- +1,, are kept at temperatures V (constant) and zero, 
respectively, for t > 0; the initial temperature of the whole being zero. As remarked in 
Sec. 1 we wish to find the quantities of heat Q, and Q, which cross unit area of the sur- 
faces x Oand zx lL, +--+ +1, up to time t. 

Let v, be the temperature in the r-th layer, and write 6, for its Laplace transform, 


that is 
b, -/ ev dt. 


/ 1/2 fe 
dr = (p a,) ’ ee J af ite ’ nN, (5) 


Then, writing 
it is found by the usual Laplace transformation precedure® that 
i, = C, sinh g,2 + D, cosh qx, 0<2z2<l, 
D, C, sinh g(a — 1, — +++ — U-,) + D, cosh q(x — +++ — b-1), (6) 


Lteee th.,<ar<ht--- +1. 


The 2n constants C, , --- , D, are to be determined from the conditions at the bound- 


aries of the layers, namely, 


vo, = V/p, x = 0, (7) 

Vv, = Dra 5 rz=1l+ih44+---+2, (8) 
bm hy Et, pe thter th, 

| | (9) 

d, = 0, r= 1l+-->- +1, 


3Cf. Carslaw and Jaeger, loc. cit., Art. 114. 
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Using the abbreviations 





&=4¢,l,, y=i],-++ ,n, (10) 
k ' Qa, 
Co. = a r=l1,--- n—Il, (11) 
ka 
the 2n equations (7) to (9) become 
dD, V/p 
C. sinh & + D, cosh &, D.., = 0, r=i,-+-,n—1| 
i (12) 
C. cosh —&, + D, sinh & g,C.4, = 0, r=I1,---,n-—1 
C,, sinh &, + D,, cosh &, 0. 


In the next section we derive some properties of the determinants which appear in 
the solution of these equations, and in Sec. 4 the solution of the present problem is 
completed. 

3. Properties of the determinants. Let A, be the determinant of 2r rows and columns 
defined by 


0 I 0 0 0 
sinh £, cosh é, 0 — | 0 ree 
cosh £, sinh &, —7O; 0 0 
0 0 sinh £ cosh é, 0 —] 
0 0 cosh & sinh &£, —3O, 0 iene (13) 
0 . 0 sinh ¢,_, cosh é,_, 0 —] 
0 . . 0 cosh —&,_, sinh &_, —o,-, 0 
0 . . 0 0 0 sinh & cosh &, 


so that A, is the determinant of the system of equations (12). Also let 2, be the determi- 
nant obtained from A, by interchanging sinh £, and cosh &, in the last row. Finally let 
6, and w, be the determinants obtained from A, and Q, , respectively, by interchanging 
0 and 1 in the first two columns of the first row. 

All the solutions in which we are interested can be expressed in terms of the four 
determinants A, , 2, , 6, , and w, . To evaluate these we proceed as follows: expanding 
A, and 2, by Laplace’s development in terms of the last two columns we find 
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A, = —A,-_,0,-, cosh &, — 9,_, sinh &, 
(14) 
2, = —A,_,¢,-, sinh — — Q,_, cosh &, , 
also A, = —simh & , 2, = —cosh &, , (15) 


and so A, and Q, can be written down by repeated application of (14). For example, 


A» a, cosh & sinh £ + sinh & cosh &, , 
2. o, sinh & sinh €, + cosh & cosh &, , 
(16) 
A, = —o 0, cosh & cosh & sinh —; — o, cosh & sinh & cosh: , 
— o, sinh é, sinh & sinh & — sinh &, cosh & cosh &, . 


Clearly, from their definitions, 5, and w, satisfy the same recurrence relations (14) as 
A, and Q, , but in place of (15) we have 
6, = cosh &, , @, = sinh é, . (17) 
To obtain the results we need, the explicit expressions for these determinants are 
not required, but only the first two terms of their expansions in ascending powers of p. 
Writing 
1/2 ‘ 
& = 7p ’ (18) 
so that by (10) 
7 = les”. r= :, PF. Ty (19) 


we find from (15) 


/2 1 
—ap'*(1 a 6 mp -}- o* ) 


A, — 

(20) 
l v4 
C = ae " ee 
2, (1 + 2 mp + ) 
We now assume 
A, = (-—1)'n,p'*(A, + Bp + --:) 

(21) 


Q, = (-1)'(A/ + Boyt :::), 
substitute these values in the recurrence relations (14), and compare coefficients. This 
gives \ 


Ae = 9-107r-14,-1 + QrAles 


34” 


n-B, = 1,-10,- 1B,-1 + ; 1:Nr—-19 rr Aya + n,-Bt_, + ; NrAr-1 
f (22) 


—_ 
~~ 


A ch 


i 
97 


247 
r4ir—i * 


Bi = 0-12 Mr—1 Ay + Bi, + 











192 J. C. JAEGER [Vol. VIII, No. 2 


This chain of equations is to. be solved successively starting with the values given 
by (20), viz. 


I : / , I 2 q 
A, =1, B, = Gm» Aj = 1, Bi =5m- (23) 
The results may be simplified a little by introducing the quantities 
g r r 
R=7, R= > R., H,=loc, Hi= dOu., (24) 
Up s=1 s =] 


which have important physical interpretations: R, and H, are the thermal resistance 
and heat capacity, respectively, of unit area of the 7-th layer, and R/ and H? are the 
corresponding quantities for the first r layers taken together. Using the notation (24) 


we get finally 


A, = R!/R. , (25) 

A; = l, (26) 

] = 

B, = (R,-1/R)B,-1 + 5 RtsH, + = Re, + Bes, (27) 

B! = R!_,H, + 5R,H, + Bi. . (28) 

(27) and (28) are to be solved successively, beginning with B, = (1/6)R,H, , Bi = 


(1/2)R,H, . The solution of (28) can be written down immediately, it is 
1 r r-l ; 
Br = 5 > R.H, + > RIAs ; (29) 


but there does not seem to be a simple expression for B, . 
In the same way, to determine the first terms of the expansions of 5, and w, we 


assume 
" (ke p c ee 
wiwtl I eae 
6, ( 1) (ky p:c,)'”* (a, + ),p + ) 
(30) 
r—1 p' ( , bh’ . 
oy = (—1) (k,p,¢,)'”* a, + ),p + ); 
and we find 
a, = 1, (31) 
at = H], (32) 


and b, and 0b} are given by 


oo. RH, 4+ RH, (33) 


bo 
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ak ek as | R.H? + 5 RH H!. (34) 


) 
with b, = 5 ul, . bi = ; RH; . (35) 
) 


It follows immediately from (33) that 
r r—-1 
b= 5 DRM. + LRH! (36) 
s=1 s=1 


4. Solution of the problem of Section 2. The solution of the equations (12) can now be 
written down in terms of the determinants of Sec. 3. We shall require C, , C, , and D, . 


It follows immediately that 


V6, 
eo (37) 


1 


Also, evaluating the determinants which occur in the numerators of C,, and D, by 
Laplace’s expansion in terms of their last two columns, we find 


” cos “sinh 
n-1 V cosh &, D, = (—1)" V sinh &, | (38) 


C, = (—T , 
) pA, pA, 


We have to find Q, , the total amount of heat which crosses unit area of the plane 
0 up to time ¢, and Q, , the total amount crossing unit area of the plane x = 1, + 


nt av, 
—k, | Bail dt, 


r= 


--+ + 1, up to time ¢. 
Q, is given by 


and its Laplace transform is 





Therefore, using (37), 





nD = kin V6, | 

Q, = - “te (39) 
Similarly, using (38), 

A = (1) ka . 

Q, = (-)) DA, (40) 


Q, and Q, are found from their Laplace transforms by the use of the inversion theorem 


followed b V contour inte gration. Thus 
Q _ 2 | e”' ( ) d R( ) > 0 
0 . 0 Pp P; Y . 


y-io 


The integrand has an infinite number of poles on the negative real axis which give 
rise to terms involving negative exponentials in the time, also it has a double pole at 
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the origin which gives the asymptote we require. To find the residue at the origin the 
expansions of 6, and A, developed in Sec. 3 are needed. Using (21) and (30) we get 


- Via, 4 bp f +++) ¥ 
Q, e’ , ” . I 
" R,p'(A, + B,p 4 eee)” (41) 


Therefore the contribution of the pole at the origin to Q, is 
vs (B )} 
Q, jt — ( - — b,)), (42 
ve = FP | A, "If 
where the values (25) and (31) of A, and a, have been used. Ff is, of course, the thermal 


resistance of the whole wall so that (42) corresponds to steady flow for time 


B,, 
l (4 — f,}. 


n 


In the same way it follows from (40) that for large values of the time 


V ( Bs) ; 
) 7 = : : 
Q,, R l re (43) 


n 


Thus the time lags at the first and last surfaces are respectively 


> b, and (44) 


These are calculated from the formulae of Sec. 3: numerical examples of the pro- 
cedure are given in See. 7. 

It remains to discuss the way in which the linear asymptotes (42) and (43) are ap- 
proached. For a complete solution it is necessary to find the zeros of A, , qua function 
of p, from which the exponential terms in the complete solution are obtained. Now 
from its definition (—A,/B,) is a crude first approximation to the smallest zero of A, , 
and, moreover, one which is certainly too small in magnitude. It follows at once that 
the departures of the curves for Q,) and Q, from their asymptotes have an exponential 


factor which decreases more rapidly than 
exp (—A,t/B,). 


In the statement of Sec. 2 there was supposed to be perfect thermal contact at the 
surfaces of the layers, and the two outside surfaces were supposed to be at temperatures 
V and zero. Normally there are air gaps or films or contact resistances of some sort at 
the outside surfaces and between the layers: these may be taken into account by re- 
garding each contact resistance as an additional layer of the appropriate thermal re- 
sistance and zero heat capacity. 

If the wall is initially at temperature V instead of zero, Q) and Q, have the same 


values as those for the wall initially at zero with x = O maintained at zero and x = 
lL, +--+. +/1,at —V fort > 0. Thus the time lags at the two surfaces are those, obtained 


as above, for flow through the wall in the reverse direction. The calculation of these is 
discussed in Sec. 6. 

Finally, it may be remarked that, since explicit expressions for 4, and 6, such as 
(16) can be written down from (14), explicit formulae for the Laplace transforms of 
(, and Q, and of the temperatures in the slab can also be written down. The solutions 


the 
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for harmonic surface temperature also follow in the usual way from these Laplace 
transforms. 

5. Other problems on the composite slab. The solutions of a number of other problems 
on the composite slab behave in the same way as Q, and Q, of Sec. 4, and in all cases 
the time lags can be expressed in terms of the eight quantities A, , --- , bf defined in 
Sec. 3. 

In all the cases below, the slab and notation are those of Sec. 2 and the slab is initially 
at zero temperature; only the boundary conditions at the surfaces are changed. 

(i) Constant flux F at x = O fort > 0. Zero temperature atx = l, + --- + 1,. 
In this case the quantity of heat Q, which has crossed unit area of the plane x = I, + 
+ /, up to time ¢ is ultimately 


Q. = Fit — b,). (45) 


(ii) Constant flux F atx = 0 fort > 0. No flow of heat atx = 1, + --- +1,. 
Here the temperature at « = 0 is ultimately 
F b/ 
7\t— G At 2 6 
fd { a, +8 (46) 
and the temperature at x = 1, + --+ + l, is ultimately 


F b; / 
/ t _— mt ° 4 
H rd (47) 
(iii) Linearly increasing temperature Vt at x = 0, and no flow of heatatx = 1,+---+1,. 


The temperature at « = 1, + --- + 1, is ultimately 
Vii — Bi} (48) 


6. Flow through the same wall in opposite directions. Most composite walls are not 
symmetrical, so that, though their thermal resistance is independent of the direction of 
flow of heat, their transient behaviour may be expected to be different for the two 
directions. It is also necessary to study flow in the reverse direction for the problem 
mentioned at the end of Sec. 4. 

Considering again the wall of Sec. 2, suppose now that the surface x = 1, + --- +1, 
is maintained at V for ¢ > 0 and the surface x = 0 at zero, the initial temperature of the 
wall being zero. The only modification is that the first and last of equations (12) are 
replaced by 

D, = 0 
(49) 
C,, sinh &, + D, cosh &, = V/p 


If Q, is the flow per unit area from right to left across the surface x = 0 up to time t, 
and Q, that across the surface x = 1, + --- + 1, , we find on solving that 


- n Venn . 
Qo = (—1" Sa (50) 


( Vk guQ, 


Q, = ah (51) 
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Since (50) is identical with (40) it follows that the quantity of heat passing through 
the wall up to time ¢ is independent of the direction of flow. This is a special case of a 


general reciprocal theorem. 
The quantities of heat passing into the wall, of course, are different, but a simple 
result about the time lags follows at once. From (51), in the usual way, we find that 


— i ( a (Be o Bi) (52) 


Comparing this with the corresponding quantity (42) it appears that b, has been 


for large values of the time 


replaced by Bj . 
These results, so far as time-lags are concerned, may be summarized by the statement 


that for flow in the reverse direction in a given wall with zero initial temperature, 
B,/A, and b{/a, are unaltered, while b, and Bf are interchanged. The latter remark 
follows also from the formulae (36) and (29). It may also be remarked here that b, is 
the time taken for a quantity of heat equal to the steady state heat content of the wall 
with constant and zero surface temperatures to flow through the wall under these condi- 
tions, while Bf is the corresponding quantity for flow in the reverse direction. 

7. Numerical calculations for a typical wall and discussion of results. First we collect 
for reference the formulae derived earlier. The fundamental thermal quantities for the 
layers are 

4 - ; . 
R=7, R= >~R., H.=loc,, Hi= >Od.. (53) 
Vr 8 l s l 
The time lags of §§4, 5 all involved the four quantities 
B./ As; B;, b, , bi/as, 
where 


A, = fe./h..« a. = Ff’ 


B.=2RH,, Be=3RH,, bw =5RH,, wf =2RHi (58) 
) “ “ ) 


>} 


B, = (R,_,/R,)B,-1 + ; R!_,H, + 2 RH, + Bi. (56) 
B! = R!_,H, + > RH, +B, (57) 
b, = b+ 3 RH, + RH. (58) 

(59) 


bf = bf_, + b,_,.H4, + : R,H, + : Ri fla: « 


We first study in detail the wall described in Table I: units are lbs., ft., hrs., B.T.U., 


and °F. 
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TABLE I, 
Layer Substance l pc k H, R, 
1 | — Air film 0 0.25 
2 Brick 0.333 25 0.42 8.33 0.79 
3 Cavity 0 0.75 
4 Brick 0.333 25 0.42 8.33 0.79 
5 Plaster 0.083 15 0.25 1.25 0.33 
6 Air film 0 0.60 
| 





For numerical calculations it is most convenient to fill in the thermal quantities H, 
and F, in a table such as Table II, then to calculate the subsidiary quantities such as 
Hi , 3R,H, , ete. from them, and finally to work down the columns for B,, B’ , b, , bf 
using (56) to (59). If only a single time lag is needed, all these columns need not be filled 




















in. 
TABLE II. 

- — — aes aime ——$—— a 
© IR 1 
,idk,. ok. Re i... Re 9 HH, A B. B b, a, bf 
—_ } | r 

, ' aa, Se eae ee | 
1 |0 0 0.25 0.25} 0 1 0 0 | 0 0 0 
2 |8.33 8.33 0.79 1.04} 0.22 2:06 0 3.29 |1.82 2.14 5.37) 3.29 8.33 9.13 
3 10 8.33 0.75 1.79} 1.05 0 6.25 0 2.39 7.62 5.37| 9.54 8.33 9.13 
4 /8.33 16.66 0.79 2.58) 0.95 14.91 6.58 3.29 |8.27 21.16 23.57/19.41 16.66 125.1 
5 |1.25 17.91 0.33 2.91] 2.39 3.22 5.50 0.21 |8.82 75.91 27.00)25.11 17.91 152.9 
6 0 17.91 0.60 3.51) 0.55 O 10.75 O 5.85 68.75 27.00/35.86 17.91 152.9 





The formulae and results for the time lags in the problem of §§2-4 and in those of 
Sec. 5 are summarized in Table III where the corresponding formulae (55) for a single 
layer of resistance R and heat capacity H are also given for subsequent reference. The 
results for the wall described in Table I are given under “‘Wall 1’’. “Wall 2” differs from 
Wall 1 only in that the layers of brick are replaced by layers of wood of thickness 0.066, 
pe = 10, k = 0.083, H, = 0.66, R, = 0.79. Wall 2 thus has the same overall thermal 
resistance R{ = 3.51 as Wall 1, but its overall heat capacity H¢ is 2.57 in place of 17.91. 

Walls 1 and 2 are idealized building walls of brick and timber construction respec- 
tively. The flow through these up to time ¢ (hours) per unit temperature difference 


between the surfaces is, for large ¢, 


_! 
3.5 


B.T.U. respectively, showing clearly the important effect of the difference in heat ca- 


(t — 1.6) 


1 
_ 8 $ — 
i (t 11.8) ind 3 Bl 


pacity of the two walls. 
This investigation was begun in connection with the flow of heat in refrigerated rail 
cars.* These are cooled by ice carried inside, and the quantity Q, of Sec. 6 determines 
*I am indebted to Mr. E. W. Hicks of the Division of Food Preservation and Transport, Council for 
Scientific and Industrial Research, Australia, for suggesting the problem and much useful discussion. 








198 J. C. JAEGER [Vol. VIII, No. 2 


TaBLe ITI. 


| | 
| | 


Boundary -_ | Single ; toes 
Problem wed Quantity | Timelag | ingle Wall 1 | Wall 2 
conditions , | layer | 
| | : a | 
| 
Sec. 2 Const. Temp. | aF | (B,/A,) — b, — ; RH | — 24.11 —2.13 
‘ 
| 
— | — 
Zero ‘Temp. Q,, B,,/A, | 6G RH | 11.75 1.59 
| : ) : 
Sec. 5(i) Const. Flux | | 
Zero ‘Temp. :. b,, | 5 RH | 35.86 | 3.42 
—- 
See. 5(ii) Const. Flux Vo | (b//a’) — BI | sate RH | —18.46} 4.11 
| | _ | | 
| a | 
Zero Flux Vp b, /a;. 6 RH 8.54 1.19 
| | ) | | 
_ —_ —————— — - 
Sec. 5(iii) | Linear Temp. | | | | 
| | 
| | | 
Zero Flux | B. 5 RH | 27.00 5.30 


the quantity of ice melted up to time ¢. The time lag in this quantity is B,, 
knowledge of this provides sufficient information to decide on the performance of a car 
in practice. A discussion of the magnitude of this time lag for various car designs will 
be given elsewhere. 

Here it is desired to stress the general applicability of these time-lags for providing 
an easily calculated criterion for comparing the performance of different walls under a 
sudden change of external conditions. Any of the cases of Table III may be used, but 
the most generally useful is B,/A, , the time lag in transfer of heat through a wall when 
its outer temperature is suddenly raised. 

Finally, the question of how far it is possible to define an 
wall’? which will have the same transient behaviour as a given composite wall should be 
discussed. If the thermal resistance and heat capacity of this equivalent homogeneous 


wall are R and H it follows from Table III that we must have 


“equivalent homogeneous 


= 2b, = 2B’ = -—; 
a, 


6B,, 6b, : 
RH = . (60) 
In actual walls, the four quantities on the right of (60) are often far from equal; thus 
for Walls 1 and 2 of Table III they are respectively 70.5, 71.7, 54, 51.2 and 9.5, 7.4, 10.6, 
7.1. Thus if an equivalent homogeneous wall is defined based on any particular transient 
phenomenon it may be expected to be a good deal out in the estimates of others. 
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—NOTES— 
AN ELECTRICAL NETWORK WITH VARYING PARAMETERS* 


By C. P. GADSDEN (Tulane University) 


2 


The aim of this paper is to study the free vibrations in a simple linear system whose 
2 parameters may vary with time. 


, R(t) PD 


a= (4.9) 
4 — C(t) ‘< q 

















Fia. 1 Fic. 2 


The system we will consider is the electrical circuit of Fig. 1. The fundamental rela- 
tions are 


ee 
@=-h=L dt? (1) 
a Ri+C'q = 0, (2) 


where ® is the flux linkages, g the condenser charge, and i the current.’ The energy 
stored in the system is E = }3C™'g’ + 4$L7'@’, and its rate of change is found to be, 
using (1) and (2), 


dE "aoe 2 
an 3° TG - 3" Sa (3) 


From (3) it is clear that if C or L is varied, energy is added to or extracted from the 
circuit. This fact has been explored by Barrow,’ Erdelyi,* Minorsky* and others to 

*Received Sept. 19, 1949. Revision received Dec. 28, 1949. This investigation was carried out under 
a contract with the Office of Naval Research. 

‘Note that the self-induced voltage is taken to be d(Li)/dt, not Ldi/dt. The former assumption 
was made by L. A. Pipes, Operational and matrix methods in linear variable networks, Phil. Mag. 25, 
585-600 (1938). The latter was made by Barrow (Ref. 2, p. 212) and by Minorsky (Ref. 4), and is 
incorrect. 

2W. L. Barrow, On the oscillations of a circuit having a periodically varying capacitance, Proc. I. R. E. 
22, 201-212 (1934). 

3A. Erdelyi, Ueber die freien Schwingungen in Kondensatorkreisen mit periodisch veraenderlicher 
Kapazitaet, Ann. Physik (5) 19, 585-622 (1934). 

4N. Minorsky, On parametric excitation, J. Franklin Inst. 240, 25-46 (1945). See bibliography. 
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show that the vibrations may be maintained, or even built up, in spite of the dissipation 
Ri’, by giving C or L a suitable periodic variation. 

It turns out to be convenient to discuss the behavior of the circuit in terms of the 
quantities g and ® together. To this end we eliminate 7 from (1) and (2) and consider 
the system of equations 


GG sn 
° hi L'®, 
(4) 
dd _ _(-1, — pr-le 
dt — ( q RL Pp, 


It is assumed that L(t), C(t) and R(t) are continuous functions of t, with continuous 
derivatives, whose values remain between some positive bounds: 0 < m < L(t), C(d), 


R(t) < M. Under these assumptions the coefficients L~’, — C™' and —RL™ in (4) 
are continuous; this is sufficient to assure that any pair q(t), ®(t,) of initial values 
uniquely determines a solution q(t), @(¢) of (4); in particular, if q(t.) = ®(t,) = 0, we 


have the trivial solution q(t) = #(t) = 0. 

Let us picture a non-trivial solution of (4) in the g, ® plane (Fig. 2). From the pre- 
ceding remark the solution curve never passes through the origin. The shape of the 
curve depicted can be justified as follows. From (4) we obtain d@/dg = —LC™'(q/®) — R. 
Thus the slope of the curve is negative in the first and third quadrants.” Next, consider 
the polar angle 8, which we will think of as changing continuously, not as reduced modulo 
2r. To be precise we define the function @ = 6(t) to be the unique continuous function 
of ¢ satisfying the relations ®= p cos 0, q = psin 6, and 0 < @(t,) < 22, where p = 
(q° + &°)'” > 0. For the rate of change of @ we have the expression 


16 ° ae sake ies ae 
7 = L' cos’ 6+ RL“ sin 6 cos 6 + C™ sin’ 0. (5) 
c 


In the first and third quadrants sin @ cos 6 > 0 so that d0/dt > M™~' cos’@ + M™ sin’ 
@ = M™ > 0. Hence the direction of increasing time is as shown by the arrow. This 
inequality shows also that the point g, ® cannot remain continually in the first or third 
quadrant longer than a time 3Mz. 

1. Oscillation. Two distinct types of behaviour are possible for a solution curve as 
t —>: (i) it may wind continually around the origin, so that q(t) and #(t) have an 
oscillatory character (about the value zero), or (ii) after a certain time it may remain in 
the second or fourth quadrant, so that g(t) and (t) become non-oscillatory. With regard 
to these two possibilities we establish several properties. 

I. Whether the vibrations are of the oscillatory or of the non-oscillatory type is inde- 
pendent of the initial values of charge and flux (provided not both are zero). This follows 
from the fact’ that the zeros of corresponding components of any two linearly inde- 
pendent solutions of (4) separate one another, while such zeros of non-trivial linearly 
dependent solutions coincide. Thus either all such solutions are oscillatory or else all 
are non-oscillatory. 


‘Note that the slope at the point g, & depends, in general, upon ¢t; thus the diagram of Fig. 2 differs 
somewhat from the usual “‘phase-plane”’ diagram. 

6M. Bocher, On the real solutions of systems of two homogeneous linear differential equations of the 
first order, Trans. Amer. Math. Soc. 3, 196-215 (1902). 








ion 


the 
der 


1950] C. P,. GADSDEN 201 


We give now an estimate of the number of complete oscillations of q(t) or (t) that 
occur in an arbitrary time interval (¢, , ¢.). This number can be measured by the quantity 
(2r)~'[A(t.) — @(t,)]. However, instead of @ we use the angle y defined by y = arc tan 
(rq/®) and the requirement that | y — @| < 432, where r = r(t) is any positive function 
of ¢ with a continuous derivative. It is clear that 6 and y always lie in the same quadrant, 
so that the quantity (27)~'[y(t.) — Y(t,)] will also serve our purpose. Similar to (5) we 
have the expression 

lr 


| d 9 . e @ 
“ = rL' cos y + (aa +r *) sin y cos ¥ + r°'C™ sin’ y 


or in terms of 2y, 


d2y 1 -17)-1 3 -17y-1 -1 -_dr\ . 
—™ rb +r oC + (rh — r'C") cos 2¥ + (RL +r di sin 2y. 


From the last relation we have the inequalities 


dy 
A(t) < GS BO, (6) 


where 
2 2°]1/2 
A(t), Bi) = : an +r'Cc' = | (v0 — rc”) + (ro +r “) | \ 
2 | dt 


Integrating (6) between ¢, and t, we have 
2ra < Y(t.) — W(t) < 2b, 
where 


a= (2x)? [ACD dt, b = 2) [BUD at. 


Hence: 

II. The number of complete oscillations of q(t) and ®(t) between the times t, and t, is 
at least a but does not exceed b. Also, if lim a = © as t, ©, the vibrations are oscillatory, 
while if lim b <@ ast, © they are non-oscillatory. 

In using the preceding result the function r(t) should be chosen to give the best 
estimate. This choice can affect the values of a and b considerably. For example, if 
L = 1,C = 1/4,R = 0, the choice r = 1 gives 1 < dy,/dt < 4, whereas if we take r = 
L'?C"'? = 2 we have dy,/dt = 2, for which the error is zero. The question of determining 
in general the best r(t) leads to a problem in the calculus of variations; however, the 
Euler equation seems to be unmanageable. 

We can of course get oscillation and non-oscillation criteria by reducing the system 


(4) to a single second order equation of the form 


TY 4 Fy =0 (7) 


ds” 


in such a way that y(s) is oscillatory or non-oscillatory if and only if g(é) or ®(é) is like- 
wise. This can be done in several ways. When it has been accomplished the classical 
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, - ° ° °,° 7, . 
comparison theorems of Sturm and other oscillation conditions’’* for (7) can be applied. 
The following result illustrates the procedure. 


Ill. If F(t) = LC™ — (R/2)’ — (L/2) dR/dt > k > O fort > T, where k is a con- 
stant, the vibrations are oscillatory, while if F(t) < 0 fort > T they are non-oscillatory for 
t > T. To show this first eliminate ® from (4) and obtain 


d aq) dq 1 : 
: ES sce Cc = 0. 8 
“( dt +h dt + 4 o (8) 
Next introduce a new independent variable s defined by s = {7 L~* dt. Since ds/dt > 0, 
t is also a single-valued function of s, say ¢ = g(s). Also, s > M~'(t — T), so that s > 
as t >, Putting g(t) = g{g(s)] = q(s), and similarly R() = R,(s), L(t) = L,(s), 
C(t) = C,(s), (8) becomes 


lq, 


dq, +R, ss 
ds 


. + L,Cy'q, = 0. (9) 


Finally in (9) put 


a8 


| 
qi(s) = y(s) exp | — 5 f(s) ds (10) 


and get an equation of the form (7) where 
f(s) = L,Cy* — (4R,)? — 4dR,/ds = LC™ — (4R)’ — (4DdR/dt = Fit). 


A comparison theorem” of Sturm applied to (7) shows that y(s) is non-oscillatory as 
s—o if f(s) < 0 for s > 0, that is, if F(t) < 0,t > T. But if y(s) is non-oscillatory, 
(10) shows that g,(s) = q(¢) is also. On the other hand, if f(s) > k > 0 this theorem shows 
that y(s) oscillates at least as rapidly as the solutions of (7) with f(s) = k. But if y(s) is 
oscillatory, so is q(t). 

IV. If the vibrations are non-oscillatory, they are also transient, that is, lim q(t) = 
lim &(t) = 0ast—~o. To prove this note first that if the vibrations are non-oscillatory, 
q(t) and ®(t) have constant and opposite signs after some time 7’, inasmuch as q, ® 
cannot remain indefinitely in the first or third quadrants of Fig. 2. Thus for ¢ > T' the 
functions | q(t) | and | ®(¢) | satisfy the system 


d|q| 


qt _ ae | ® |, (11a) 
lel 
d ia =C"|q|—RL"| 4}. (11b) 


From (lla) | g(t) | decreases monotonically as ¢ increases and so approaches a definite . 
limit | gi) | ast >. From (11b) we have 


t t 
: | ®@ | exp [ RL at| = C'|q| exp [ RL” dt. (12) 
dt L Jr JT 

7R. Bellman, A survey of the theory of boundedness, stability and asymptotic behavior of solutions of 
linear and non-linear differential and difference equations, Office of Naval Research, Department of the 
Navy, Washington, D. C., 1949, p. 116-136. 

8A. Wintner, A criterion of oscillatory stability, Q. Appl. Math. 7, 115-117 (1949). 
°F. L. Ince, Ordinary differential equations, Dover Publications, New York, p. 226. 


ad. 


N= 


or 
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Integrating (12) between 7 and ¢ gives 


| &(¢) | = [ #(7) | + l C"'lq| (exp [ RL” at) at|(esp | RL" at) : (13) 
“T T 4 


To the right side of (13) we apply an extension’® of L’Hospital’s rule concerning the 
limits of the quotient f(t)g"*(t) as t >: if lim | g(f)| = as t—~, and f(é) and g(t) 
are never simultaneously zero or infinite in the interval (7, ~) then 

f(t) f' 


lim sup —~ < lim sup 3 
oP g(t) "ve 





and 





——s ao» Fie 
lim inf a(t) 2 lim inf q(t)’ 


where the primes denote derivatives. In our case 


at 
g(t) | = exp | RL™ dt > exp mM"“'(t — T) as to, 
- 


o 


Hence we have 


-1 
lim sup | &(t) | < lim sup os < lim sup m™*M | q(b) | 


(14) 
= |q(~)|m°*M <a, 
and 
-1 
co> lim inf | &(é) | > lim inf ote] > lim inf M~*m | a(é) | 
(15) 


= |q(o) | Mm >0, 


where all limits are taken as i ~«. Now suppose lim inf | 6(¢) | > 0. Then from the 
definition of the lim inf, if « is a number such that 0 < « < lim inf | ®(é) |, there is a t, 
such that for ¢ > t, , | ®(t)| > «. Hence from (lla) d | q|/dt < —M~‘e. Integration of 
this inequality gives 0 < | g(t) | < | g(t) | — M™'e(t — 4). But as t > the right side 
of this inequality — — . This contradiction shows that lim inf | (é) | = 0. Then from 
(15), | g@(@) | = 0 and from (14), lim sup | &(é) | = 0, so that lim q(é) = lim #() = 0 
ast—o., 

2. Stability. In this section we prove some facts concerning the stability of the 
solutions g(t), &(#) of (4). It should be emphasized that unstable solutions are indeed 
possible. This has been shown both theoretically and experimentally.?"*"* 

We will say that a function h(t) is stable if it remains bounded as t >, i.e., if lim 
sup | A(t) | <. Our first result shows that g(¢) and #(t) have the same stability be- 
havior. The proof requires the assumption that dR/dt remains bounded as t —@. 


1o],, M. Graves, The theory of functions of real variables, McGraw-Hill Book Co., Inc., New York, 
1946, p. 83. 
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comparison theorems of Sturm and other oscillation conditions’’* for (7) can be applied. 
The following result illustrates the procedure. 


III. Jf F(t) = LC™* — (R/2)’ — (L/2) dR/dt > k > O fort > T, where k is a con- 
stant, the vibrations are oscillatory, while if F(t) < 0 fort > T they are non-oscillatory for 
t > T. To show this first eliminate ® from (4) and obtain 


d aq) dq a 
= i= —* + C = (0. 8 
dt (x dt +h dt + q ° (8) 
Next introduce a new independent variable s defined by s = f7 L™ dt. Since ds/dt > 0, 
t is also a single-valued function of s, say t = g(s). Also, s > M~'(t — T), so that s ~@ 
as t >. Putting g(t) = g{g(s)] = g(s), and similarly R() = R,(s), Lit) = L,(s), 
C(t) = C,(s), (8) becomes 


4h +p, + 1.07'q, = 0. (9) 
ds ds 
Finally in (9) put 
om [ Tis 
qi(s) = y(s) exp — 5 Rs) ds (10) 


and get an equation of the form (7) where 
f(s) = L,Cy" — (3R,)’ — 4dR,/ds = LC™ — ($R)’ — (4D)dR/dt = F(t). 


A comparison theorem’ of Sturm applied to (7) shows that y(s) is non-oscillatory as 
s—o if f(s) < Ofors > 0, that is, if F(t) < 0,¢ > T. But if y(s) is non-oscillatory, 
(10) shows that ¢,(s) = q(é) is also. On the other hand, if f(s) > k > 0 this theorem shows 
that y(s) oscillates at least as rapidly as the solutions of (7) with f(s) = k. But if y(s) is 
oscillatory, so is q(t). 

IV. If the vibrations are non-oscillatory, they are also transient, that is, lim q(t) = 
lim @(t) = 0ast—~o. To prove this note first that if the vibrations are non-oscillatory, 
q(t) and ®(t) have constant and opposite signs after some time 7’, inasmuch as g, ® 
cannot remain indefinitely in the first or third quadrants of Fig. 2. Thus for ¢ > T the 
functions | g(t) | and | ®(¢) | satisfy the system 


diq| _ _ 7-1; 9) ‘ 
: lian LL | |, (11a) 

d ® | iy »7T-!1 

~_ =x ( q|— RL | @|. (11b) 


From (lla) | g(¢) | decreases monotonically as ¢ increases and so approaches a definite . 
limit | g(@) | ast ~o. From (11b) we have 


t t 
_ | ® | exp [ RL” at| = C'|q| exp [ RL” dt. (12) 
JT JT 
7™R. Bellman, A survey of the theory of boundedness, stability and asymptotic behavior of solutions of 
linear and non-linear differential and difference equations, Office of Naval Research, Department of the 
Navy, Washington, D. C., 1949, p. 116-136. 
8A. Wintner, A criterion of oscillatory stability, Q. Appl. Math. 7, 115-117 (1949). 
°E. L. Ince, Ordinary differential equations, Dover Publications, New York, p. 226. 
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Integrating (12) between 7 and ¢ gives 


| &(t) | = [ a(7)|+ [ c*1a| (exp [ RL” at) at|(exp [ RL” a) (13) 
oF T ig 


To the right side of (13) we apply an extension” of L’Hospital’s rule concerning the 
limits of the quotient f()g"*(t) as t >: if lim | g(t) | = as t—o, and f(t) and g(t) 
are never simultaneously zero or infinite in the interval (7, ~) then 

I) f' 


lim sup —~ < lim sup 
"a os ee 





and 


Lim inf i > lim inf : om 


where the primes denote derivatives. In our case 





g(t) | = exp [ RL™ dt > exp mM“'(t — T) — © as to, 
JT 
Hence we have 


-1 
lim sup | &(d) | < lim sup cis | < lim sup mM | a(é) | 


aL (14) 
= | q(~) | m’°M <o, 
and 
-1 
co > lim inf | &(t) | > lim inf la | > lim inf M~’m | q(é) | 
(15) 


| q(@) | M*m > 0, 


where all limits are taken as t ~ «©. Now suppose lim inf | @(¢) | > 0. Then from the 
definition of the lim inf, if « is a number such that 0 < « < lim inf | (é) |, there is a ¢, 
such that for t > t, , | &(t) | > «. Hence from (lla) d | q|/dt < —M™*e. Integration of 
this inequality gives 0 < | g(t) | < | g(t,) | — M~'e(t — t,). But as t © the right side 
of this inequality — — . This contradiction shows that lim inf | @(é) | = 0. Then from 
(15), | g(~) | = 0 and from (14), lim sup | (é) | = 0, so that lim g(t) = lim &(t) = 0 
ast—o., 

2. Stability. In this section we prove some facts concerning the stability of the 
solutions g(t), &(t) of (4). It should be emphasized that unstable solutions are indeed 
possible. This has been shown both theoretically and experimentally.”'*"* 

We will say that a function h(t) is stable if it remains bounded as t >, i.e., if lim 
sup | A(t) | <. Our first result shows that g(t) and #(t) have the same stability be- 
havior. The proof requires the assumption that dR/dt remains bounded as t >. 


],, M. Graves, The theory of functions of real variables, McGraw-Hill Book Co., Inc., New York, 
1946, p. 83. 
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V. If q(t) is stable or transient, P(t) is likewise, and conversely. From (4) we get 


an equation like (12) 


d at t 
| # exp | RL at | = —C™'q exp [ RL” dt, 
dt Jr 


7T 


and like (13) 


at at t -1 
(7) = | er) _ | ay ‘(exp | RL” at) at|(cxp [ RL™ at) . 
JT JT J7 


7 


at at t =] 
Clq! (exp | RL™ a) at|(cxp [ RL" iu) ; (16) 
“T es 


7 


Hence 


| (2) ra] a(T) | + | 


Applying the extended L’Hospital rule to (16) gives 
le q(t) | 


' , ( M ,. 
lim sup | &(¢) | < lim sup < —3 lim sup | q(d) |. 
m 


PL 


Thus if g(t) is stable, so is &(t), while if g(t) — 0 as t ~@, so does #(t). To show the 
converse we obtain from (4) 


dq 


int, ap OO 
ath ( 7 = if 


dt 


This equation is satisfied by 


at Im nt at -1 
q(t) = | a(n) — | R e (exp | oS at) at |(cxp I R'C at) - (17) 
7 /7 


“7 


An integration by parts puts (17) in the form 


q(t) = —R'&(t) + | ar) 
pt ; ] > 4 at me - at = =-1 
+ I OR ‘(c ‘— ah \ exp ‘ ie hy it) at|(cxp | RVC at) : 


Let k be an upper bound for R~* | C~' — dR/dt| in the interval t > 7. Then 
PI 


| q(t) | < R™* | Hd) +[ la 


nt at ‘ ' 
k | (ex | ie aie u) u|( 48) i i) . 


Hence lim sup | g(t) | < (m™* + kM”) lim sup | #(é) |. Thus if &(¢) is stable or transient, 
q(t) is likewise. 

Next let us see if we have an analog for stability of the result J for oscillation. That 
is, does the behavior of the circuit with regard to stability depend upon the initial 
conditions? Strictly speaking it does. Thus one (non-trivial) solution of the system (4) 


may be stable and another unstable. In such a case, however, we can see that the system 





vet 














1950] Cc. P. GADSDEN 205 


is essentially unstable. That is, the possible initial values qo , ®) which give rise to stable 
solutions all lie on a certain straight line through the origin of the g, ® plane. For suppose 
g:(t), ,(t) is a (non-trivial) stable solution and q,(t), &.(¢) an unstable solution. These 
are clearly linearly independent, so that every solution of (4) has the form 


q(t) = & q(t) + b q2(t), 
(18) 


P(t) = a #,(t) + b (0), 


where a and 6 are constants. If the solution (18) is to be stable we must certainly have 
b = 0. Putting ¢ = ¢, in (18) and eliminating a results in 


Q(t) Po — ®, (to) Go = 0, 


a straight line in qo , Po. 
As stability criteria we have, besides IV, 
VI. All solutions of (4) are stable if either 


. d{C —— 
(i) i (2) < 2RCL 


= d{C 
(ii) #(¢) 20 


holds for t > T for some T. To prove that (i) is sufficient to insure stability consider 
the positive form E, = CL~'&’* + q’. If q, ® is a solution of (4) we have 


dE, _ seis 4 (¢) | . 
dt | anc L dt \L ms 
so that EF, is a non-increasing function of t for t > T if (i) holds. Being positive, FE, is 


therefore bounded as t >. Finally E, > mM~'#* + q’ > 0, so that q(t) and ®(t) are 
stable. To prove condition (ii) consider Z, = 6° + LC~'g’. In this case we have 


a = _| ae ‘ (5) Jo — 2RL™'®’. 
Thus if (ii) holds, Z, is non-increasing and q(t) and #(¢) are stable. The conditions (i) 
and (ii) show that the derivative of C/L must vary continually outside both the limits 
0 and 2RCL ~~ if the vibrations are unstable. 

3. Remarks. The above results apply at once of course to the dual of the circuit of 
Fig. 1, namely a parallel connection of one of each of the three types of elements, as 
well as to mechanical or acoustical analogs of these circuits. The question of the de- 
scription of the forced vibrations may be raised. As is well known from the theory of 
systems of linear differential equations, the forced component can be expressed in terms 
of the free vibrations. Thus for any discussion of the forced vibrations it should be helpful 
to know something of the character of the (relatively) simpler free vibrations. 


or 
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A METHOD FOR OBTAINING THE CHARACTERISTIC EQUATION OF A 
MATRIX AND COMPUTING THE ASSOCIATED MODAL COLUMNS* 


By HENRY E. FETTIS (Air Materiel Command) 


The method described below is a process by which the characteristic equation of a 
matrix may be established without resorting to direct expansion of the coefficients by 
minors. Further, once the roots of the characteristic equation have been found, the 
corresponding eigenvectors can be found directly without the additional labor of solving 
a set of simultaneous equations. Since the only operations involved in the process are 
the standard ones of matrix multiplication and addition, the work can be set down and 
performed in routine fashion. The process may also be adapted to punch-card methods 
for matrices of large order. 

1. Preliminary considerations. Consider the matrix 








: ‘ 
A Ay, Qy> se Ay, 
As, Aro Qo, 
(1) 
LA Gs Ann 4 
The characteristic equation of A is defined as 
ai1,—A Ayo tee Qin 
ies Aso — XA a 
== 0. (2) 
An Ane cee Ann — XA 








When written out in full, (2) is a polynomial equation of degree n in X: 
nN" — BA" + DA? + eee + OO, = O. (3) 


The » roots of (3) are called the characteristic numbers or eigenvalues of the matrix A; 
these numbers are the values of \ for which the homogeneous system of equations 


(ay, — A)X, + (Ay2)t%2 + +++ (a,,)z, = O, 


(@2,)4, + (dog — A). + +++ (Ay,)zX, = O, 
(4) 


(Gn1)%1 + (An2)X2 + +++ (Quan — ADH, = O, 


(or, in matrix form, Ax = Az) possesses a non-trivial solution. 





*Received Nov. 8, 1949. 
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The expansion of (2) in the form (3) requires the evaluation of sums of determinants 
of successively higher order, commencing with b, as the sum of the n elements on the 
principal diagonal. b, is called the trace of the matrix A, written 


b, = TrA. (5) 


Similarly, b, is the sum of all determinants of order 2 whose diagonal’ elements coincide 
with the diagonal elements of A. In all there will be n(n — 1)/2 such determinants. In 
a like manner, }, is made up of the sum of the n(n — 1) (n — 2)/6 determinants of order 
3 which can be formed in this way; continuing, it is found that b, = | A |. It is evident 
that as the order of the matrix increases, the labor involved in expanding the char- 
acteristic equation also increases, but much more rapidly, so that a point is soon reached 
for which direct expansion is impractical. 

In addition to the work required in deriving the characteristic equation, the com- 
plete solution of the problem will in general require the resubstitution of the roots of 
the equation into the original system of equations to obtain the corresponding relation 
between the x, . Since for \ equal to any root of (3), the system (4) becomes consistent, 
it suffices to solve for (n — 1) of the zx, in terms of any arbitrary one. The relation be- 
tween the x; which is now determined apart from an arbitrary numerical factor (whose 
actual value is usually of no concern) is called the modal column or eigenvector corre- 
sponding to the characteristic number. Obviously for a matrix of order 4 or higher such 
a solution would become extremely laborious in view of the fact that for each of the n 
characteristic numbers, a system of (n — 1) simultaneous equations must be solved. 

In the following section a method is set forth by which the characteristic equation 
of a matrix may be established in a routine manner employing the standard operations 
of matrix multiplication and addition. The labor does not increase with the order of 
the matrix to the extent that direct expansion does. In addition the results so obtained 
eliminate the necessity of solving a set of simultaneous equations for each root of the 
characteristic equation. 

2. Summary of the method. Let the given matrix be A, and define the successive 
matrices A, and numbers b, as follows: 


A,=A : b, =(Tr Ao) ; 
A, = b, Ao — Ao b, = (Tr A,)/2; (6) 
A,=b,Ag— AoA ,; b, = (Tr A,)/3; 
in general 
A, = b, Ay — Ao Ax-1 ; bes, = (Tr A,)/(k + 1). 


Continuing the above, it will be found that A, = 0. This result serves as a check on 
the correctness of the operations. When this point has been reached, the characteristic 


equation is 
x” — BA" + BAT? +--+ + (—-)"D, = 0. (7) 


1The term “diagonal elements” is used to denote the elements of the principal diagonal. 
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Further, if \; is any root of (7), the corresponding eigenvector x“ is proportional to 
any column of the matrix 


Addi a A,y;' + A.” + ith (—)"" Ay_aV¢ ° (8) 


In actual practice it is not necessary to compute all the elements of the various 
matrices. In fact, if n > 4, only (n — 3) of the matrices A, need be computed in their 
entirety. F 

The method and its modifications are illustrated in the next section. 

3. Numerical example’. For the matrix 

















A =[ -2 —2 0 3 —1] 
—2 0 —3 5 0 
0 —3 —5 i 1 
3 5 1 -3 -1 
od 0 I ae: | 
we have 
i,=[ -2 —2 0 30Ci‘ &é-+*d —2 
—2 0 —3 5 0 0 
0 —3 —§ l 1 =~§ 
3 5 1 —3 -1 5 
—1 0 l —1 —11 —2 
i,=[ 18 19 8 —24 0 21 
19 38 200 -24 -6 47 
8 20 36 -24 —7 33 
—-24 —-24 —24 45 2 —25 
0 —6 —7 2 41 —7 
Thus, 6, = —11. 


(The columns on the right consist of the sums of the elements of the respective rows. 
These serve as a check on the matrix multiplication, since the product of the sum column 


*This example is taken from W. Kincaid’s paper Numerical methods for finding the characteristic 
roots and vectors of matrices, Q. Appl. Math, 5, 320-345 (1947). 








to 
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of any matrix into any other matrix should equal the sum-column of the product matrix.) 




















Continuing we compute A, = b,Ay — Aj and AoA; : 
A, = A,A; = 
rf 4 3 -8 -9 I) 14; [ -—52 -—29 33 35 —14] —27 
3-38 13 —3]1 6 |—47 —29 —200 24 -81 35]-—251 
—8 13 19 138 -—4 33 33 24 —125 25 18] —25 
—9 —3] 13 —12 9 }—30 35 —8l 25 —142 25/-—138 
111 6 -4 9 7 29 L-14 35 18 25-311 33 
Thus b, = —20/2 = —10. Next we find A, from A, = b,Ay — AoA, , and the diagonal 


elements only of the product A,A2 . 














A, = Ao: As = 
[ 72 49 —33 —65 24] 47]; | — 461 
49 200 6 31 —35] 251 39 
—33 6 175 —35 —28] 85 —956 
—65 31 —35 172 —15] 88 —576 
24-35 26-15 414-13 —78 | 





Thus, b, = 660/3 = 220. The diagonal elements of A; and the coefficient b, may now be 
found from A; = };A)9 — AoA2 and b, = Tr A; . The expression for A, is rewritten as 


A, = by Ao —(bs AG — AG A2). 
Again only the diagonal elements are needed. 
A, = [21 . A’, A, = [3523 
—39 8117 
— 144 7192 


— 84 9366 














— 142 L 540 - 


Thus, ), = —388/4 = —97. 
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A,:A; = b, Ai — A}: A, = [433 1; 
243 
728 
534 
) 340! 
A, = b, Ay — A, A; = — 243 
— 243 
— 243 
— 243 
— 243 
Hence b, = —243 and the characteristic equation is 
d* + 114° — 10d*° + 220d? — 97A + 243 = 0. 
The roots® of this equation are 
A, = —9.886487, dA. = —4.75775, As = 4.22365, A, = —1.43300, A, = .85355. 


tewriting the first columns of Ay , A; , Az and computing the first columns of A, and 
A, by column multiplication gives 


(A.)=[-2]; (A)=[ 4]; (Ad)=[ 72]; (As)=b3(A0)—(A0)(A2)= [21]; 
—2 3 49 —70 
0 —8 —33 23 
3 —9 —65 61 

\—11 ate | 98) L— 156 | 


























where (A) denotes the first column of A. Finally 


8A convenient method for solving equations of this type is described by Shih-Nge Lin, A method of 
successive approximations of evaluating the real and complex roots of cubic and higher equations, J. Math. 
Phys., 20, 231-242 (1941). 





bo 





Thus for A, 


1950] 


(A,) = b,4(Ao) — Ao(As) = 
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— 243]. 
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As a further check, it is noted that the elements of (A,) not previously computed are 


zero. 


The eigenvectors are therefore (after removing a factor of \,;) proportional to 


. 





97 


bo 





| 4 


as —f 473+ 7 72) 
3 49 

—S8 —33 

~9 —65 

11 4 24 4 





19.772974 
15.772974 


19.772974 


-16.772974 


8 


0 
8 


9 


29.65946. 
20.65946 


1] 


9.886487 
1.113513 











x — | 21 | 


61 








peed 156 4 





A, + [ —243 ] 








Q - 


























+ 72 - 21 — 243 
— 155.9393 + 829.8648 — 7996.83 
— 83.9393 + 808.8648 |— 8239.83 | 
+ 49 + 70 + 0 

— 165.8258 + 1154.997 _ — 12110.917 
— 116.8258 + 1224.997  |— 12110.917 
- = —- 2 + 0 
— 79.0919 + 1108.195 — 10728.766 
— 112.0919 + 1085.195  |— 10728.766 
— 65 — 61 + 0 

+ 204.2495 — 1376.688_ _+_:14213.68 
+ 139.2495 — 1437.688  |+ 14213.68 
+ 24 + 156 + 0 

+ 11.00873 — 346.11335 + 1879.553 
+ 35.00873 — 190.11335 |4+ 1879.553 | 


A convenient method of evaluating the above polynomials is by synthetic division. 
— 9.886487, 





9 
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Therefore the first eigenvector is proportional to 


| —8239.83 | 
— 12110.92 
— 10728.77 


14213.68 








1879.55. 


) 000000 ] , 


469802 


302062 


.724997 








L— 2281064 
The remaining eigenvectors, found in a similar manner are: 


a” =[ 1.00000]; 2” =[ .04830]; 2° =[ .4224 ]; 2” =[ 1.0000) 
08337 1.00000 — .00483 — .4364 
— 1.0793 — .2703 3993 .1825 


— .6869 7018 4099 4350 


























5302 L— 1935 | 1.0000 J L— .6757 | 


In this example, direct expansion of the characteristic equation would involve the 
computation of 10 determinants of order 2, 10 of order 3, 5 of order 4, and 1 of order 5. 
In addition the only available check would be an independent evaluation. Further, for 
each characteristic number, the solution of a system of four equations would be required, 
necessitating the additional evaluation of 5 determinants of order 4 for each X,; . 




















